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1  Introduction 


General 

The  degree  to  which  aquatic  macrophytes  influence  the  ecosystem  is 
proportional  to  plant  mass  and  depends  on  plant  species  and  physicochemical 
factors.  Therefore,  predictions  of  the  environmental  impact  of  management 
measures  concerning  aquatic  communities  should  be  based  on  accurate  estimates 
of  (a)  plant  species  and  mass  and  its  pertinent  physiological  properties,  (b)  the 
contribution  of  plants  to  the  various  food  chains,  and  (c)  the  contribution  of  the 
decay  of  plants  to  biogeochemical  cycling  and  oxygen  regime.  A  simulation 
model  for  metabolism  and  growth  of  aquatic  community  types  may  serve  as  a 
useful  tool  in  this  respect. 

Although  the  number  of  simulation  models  for  growth  of  monotypic, 
submersed  macrophyte  communities  is  increasing  (e.g.,  Titus  et  al.  1975;  Best 
1981;  Collins  and  Wlosinski  1985;  Best  and  Jacobs  1990;  Hootsmans  1991, 
1994;  Scheffer,  Bakema,  and  Wortelboer  1993;  Best  and  Boyd  1996),  it  is  still 
relatively  low  compared  with  that  for  terrestrial  vegetation.  The  current  model 
has  been  developed  because  none  of  the  existing  models  were  suitable  to 
simulate  the  behavior  of  a  monotypic  milfoil  community  under  various 
environmental  and  climatological  conditions  over  a  period  ranging  from  season 
to  several  years. 


Distribution  of  Eurasian  Watermilfoil 
within  the  United  States 

The  submersed,  rooted  aquatic  macrophyte  Myriophyllum  spicatum  L.  or 
Eurasian  watermilfoil  belongs  to  the  dicotyledonous  family  Haloragaceae.  It  has 
the  ability  to  survive  unfavorable  environmental  conditions  and  has  been 
demonstrated  to  outcompete  many  other  submersed  aquatic  plant  species  in 
temperate,  subtropical,  and  tropical  areas.  This  species  has  consequently  a  very 
large  distributional  area.  It  may  be  considered  as  the  most  aggressive  member  of 
a  circumboreal  complex  of  closely  related  taxa  (Patten  1954).  A  problem  in 
discussing  the  distribution  and  rapid  spread  of  Eurasian  watermilfoil  is  that  this 
plant  species  is  morphologically  very  similar  to  the  native  North  American 
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milfoil  variously  named  Myriophyllum  exalbescens  Fern.,  M.  spicatum  var. 
exalbescens  (Fern.),  and  M.  spicatum  suhsp.exalbescens  (Fern.)  Hult.  The 
taxonomic  distinction  probably  has  not  been  made  in  all  cases  when  these  two 
species  have  been  discussed  in  literature.  Hereafter,  Eurasian  watermilfoil  will 
be  referred  to  simply  as  milfoU. 

Milfoil  is  a  native  of  Eurasia.  It  has  been  present  in  the  United  States  since 
1948  (Couch  and  Nelson  1985).  This  species  was  not  considered  a  weed  until 
the  late  1950s.  Since  that  time,  it  has  spread  from  the  east  to  the  west  coast  in 
both  the  United  States  and  Canada  (Reed  1977;  Aiken,  Newroth,  and  Wile 
1979),  and  it  has  been  documented  in  44  of  the  States  and  3  Canadian  provinces 
(Engel  1993).  Spreading  of  species  over  large  distances  was  partly  related  to 
aquarium  and  aquatic  nursery  trade  (Reed  1977).  Short-distance  dispersal 
probably  occurred  by  transport  of  plant  fragments  between  lakes  on  boats  or 
trailers  (Scales  and  Bryan  1979).  The  explosive  growth  appears  to  follow  major 
environmental  disruptions  (Nichols  and  Shaw  1986).  For  example,  the 
Chesapeake  Bay  population  increased  only  in  the  1950s  and  early  1960s  (Allen 
1973;  Bayley  et  al.  1978)  after  hurricanes  hit  the  area  repeatedly  causing 
temporarily  increased  salinity,  sedimentation,  and  inflow  of  nonpoint  source 
pollutants.  Increased  milfoil  growth  in  Cayuga  Lake,  New  York,  and  Lake 
Mendota  and  Lilly  Lake,  Wisconsin,  is  attributed  to  major  natural  or  human 
caused  disruption  (Lind  and  Cottam  1969;  Oglesby  et  al.  1976;  Nichols  1984). 
Dramatic  population  fluctuations  appear  to  be  characteristic,  since  they  have 
been  reported  not  only  in  the  native  Eurasian  range  of  milfoil  (Lundegardh- 
Ericson  1972;  Jeschke  and  Muther  1978)  but  also  in  the  Chesapeake  Bay  area 
and  in  Lake  Wingra.  In  the  Chesapeake  Bay  area,  milfoil  declined  first  in  the 
most  recently  colonized  areas  rather  than  in  the  original  epicenters  of  growth 
(Bayley  et  al.  1978),  as  such  suggesting  a  pattern  of  spreading  from  optimal 
growth  areas  to  less  optimal  ones  (Nichols  and  Shaw  1986).  Causes  of  declines 
are  still  under  discussion,  but  initial  stages  of  declines  are  commonly  attributed 
to  a  large  decrease  in  water  transparency  as  a  consequence  of  increases  in  total 
suspended-solids  concentrations  and  in  algal  growth,  respectively. 

Milfoil  is  considered  a  nuisance  plant  in  parts  of  the  United  States,  since  it 
may  interfere  with  human  utilization  of  fi'eshwater  resources,  become 
aesthetically  displeasing,  or  displace  desirable  indigenous  communities.  From  a 
shoreline  perspective,  the  biomass  in  a  dense  “mat”  of  submersed  weeds  appears 
to  be  enormous.  However,  data  on  total  biomass  and  productivity  indicate  that 
they  are  small  compared  with  those  of  several  terrestrial  plant  communities 
(Spencer  and  Bowes  1990).  This  apparent  anomaly  may  be  largely  due  to  the 
uneven  distribution  of  biomass  over  the  water  column,  with  typically 
>60  percent  concentrated  in  the  upper-water  layers. 

The  simulation  model  developed  in  this  study  concerns  Eurasian  water- 
milfoil.  The  following  appendixes  are  included  in  this  report:  Model  Listing  as 
Appendix  A,  Variable  Listing  as  Appendix  B,  and  Manipulation  of  Literature 
Data  Used  for  the  Model  Equations  as  Appendix  C.  A  user  manual  is  published 
separately  (Best  and  Boyd,  in  preparation). 
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2  MILFO:  Description  of 
Model 


Modeling  Concepts 

The  MILFO  (Version  1.0)  model  simulates  growth  of  a  typically  monoecious 
Eurasian  watermilfoil  community.  In  the  model,  growth  is  considered  the  plant 
dry  matter  accumulation  including  rhizome/root  crown  formation,  under  ample 
supply  of  nitrogen  and  phosphorus,  in  a  pest-,  disease-,  and  competitor-free 
environment  under  the  prevailing  weather  conditions.  Two  or  three  plant 
cohorts  in,  respectively,  temperate  or  tropical  areas  wax  and  wane  per  season 
with  one  and  the  same  rhizome/root  crown  system  as  a  common  basis.  The  rate 
of  dry  matter  accumulation  is  a  function  of  irradiance,  temperature,  CO2 
availability,  and  plant  characteristics.  The  rate  of  CO  2  assimilation 
(photosynthesis)  of  the  plant  community  depends  on  the  radiant  energy  absorbed 
by  the  canopy,  which  is  a  function  of  incoming  radiation,  reflection  at  the  water 
surface  and  attenuation  by  the  water  column,  attenuation  by  the  plant  material, 
and  leaf  area  of  the  community.  From  the  absorbed  radiation,  the  photosynthetic 
characteristics  of  individual  shoot  tips,  and  the  pH-determined  CO2  availability, 
the  daily  rate  of  gross  CO2  assimilation  of  the  community  is  calculated.  These 
calculations  are  executed  in  a  set  of  subroutines  added  to  the  model. 

Part  of  the  carbohydrates  produced  is  used  to  maintain  the  existing  biomass. 
The  remaining  carbohydrates  are  converted  into  structural  dry  matter  (plant 
organs).  In  the  process  of  conversion,  part  of  the  weight  is  lost  in  respiration. 
The  dry  matter  produced  is  partitioned  among  the  various  plant  organs  using 
partitioning  factors  defined  as  a  function  of  the  phenological  cycle  of  the 
community.  The  dry  weights  of  the  plant  organs  are  obtained  by  integration  of 
their  growth  rates  over  time.  The  plant  winters  through  a  system  composed  of 
root  crowns  attached  to  a  rhizome  system  in  the  sediment  with  or  without 
aboveground  plant  biomass  present.  All  calculations  are  performed  on  a  m^ 
basis.  Since  environmental  factors  and  plant  growth  characteristics  vary  with 
depth,  in  the  model  the  water  column  and  associated  growth-related  processes 
have  been  partitioned  in  0.10-m  depth  classes  (Titus  et  al.  1975). 
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Seed  formation  has  not  been  included  in  the  model  because  its  role  in 
maintaining  an  existing  milfoil  community  at  the  same  location  is  minimal 
(Hartleb,  Madsen,  and  Boylen  1993).  Dispersal  and  colonization  of  new  habitats 
by  plant  fragments  and  seeds  are  recognized,  important  characteristics  of 
Eurasian  watermilfoil.  The  latter  processes,  however,  are  better  described  using 
other  modeling  approaches  (based  on  logistic  regression  or  on  descriptions  of 
population  dynamics  varying  in  time  and  in  space),  as  discussed  by  Scheffer 
(1991). 

MILFO  requires  as  input  physiological  properties  of  the  plant  community  (in 
this  case  of  milfoil)  and  of  the  actual  environmental  and  weather  conditions  at 
the  site,  characterized  by  geographical  latitude  and  longitude,  i.e.,  water 
temperatures  (optional),  alkalinity,  pH,  and  daily  maximum  and  minimum 
temperatures  and  irradiance  for  each  day  of  the  year.  It  can  be  run  for  periods  of 
1  to  5  years. 


Modeling  Approach 

MILFO  is  a  mechanistic  model  that  explains  plant  growth  on  the  basis  of  the 
underlying  processes,  such  as  CO2  assimilation  and  respiration,  as  influenced  by 
environmental  conditions.  This  type  of  model  follows  .the  state-variable 
approach  in  that  it  is  based  on  the  assumption  that  the  state  of  each  system  can 
be  quantified  at  any  moment  and  that  changes  in  the  state  can  be  described  by 
mathematical  equations.  In  this  type  of  model,  state,  rate,  and  driving  variables 
are  distinguished.  State  variables  are  quantities  such  as  biomass  and  number  of 
individuals  of  a  population.  Driving  variables  characterize  the  effect  of 
environment  on  the  system  at  its  boundaries,  such  as  climate  and  food  supply. 
Each  state  variable  is  associated  with  rate  variables  that  characterize  its  rate  of 
change  at  a  certain  instant,  as  a  result  of  specific  processes.  These  variables 
represent  flows  of  material  between  state  variables,  the  values  of  which  are 
calculated  from  the  state  and  driving  variables  according  to  knowledge  of  the 
physical,  chemical,  and  biological  processes  involved.  After  calculating  the 
values  of  all  rate  variables,  they  are  then  used  to  calculate  the  state  variables 
according  to  the  scheme:  state  variable  at  time  t  +  At  equals  state  variable  at  time 
t  plus  the  rate  at  time  t  multiplied  by  At.  This  procedure,  called  numerical 
integration,  gives  the  new  values  of  the  state  variables,  from  which  the 
calculation  of  rate  variables  is  repeated.  To  avoid  instabilities,  the  time  interval 
At  must  be  small  enough  so  that  the  rates  do  not  change  materially  within  this 
period.  This  is  generally  the  case  when  the  time  interval  of  integration  is  smaller 
than  one-tenth  of  the  “time  coefficient”  or  “response  time.”  This  characteristic 
time  of  a  system  is  equal  to  the  inverse  of  the  most  rapid  relative  rate  of  change 
of  one  of  its  state  variables.  The  smaller  the  time  coefficient,  the  smaller  the 
time  interval  of  integration  (Rabbinge  and  De  Wit  1989). 

The  predictive  ability  of  mechanistic  models  does  not  always  live  up  to 
expectations.  It  should  be  realized,  however,  that  each  parameter  estimate  and 
process  formulation  has  its  own  uncertainty,  and  that  uncertainties  in  parameter 
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estimates  may  accumulate  in  the  prediction  of  the  final  yield.  The  primary  aim 
of  this  model  is  to  increase  insight  in  the  system  studied  by  quantitatively 
integrating  the  current  knowledge  in  a  dynamic  simulation  model.  By  studying 
the  behavior  of  such  a  model,  better  insight  in  the  real  system  is  gained. 


Implementation 

The  MILFO  model  was  implemented  as  a  FORTRAN??  program.  For 
numerical  integration,  the  Runge-Kutta  technique  is  used,  which  allows 
employing  a  variable  time-step.  The  program,  as  it  is  being  run,  integrates  the 
equations  once  per  day  in  the  main  subroutines  (MODEL,  CHRT2,  CHRT3;  see 
Figure  1),  once  per  second  in  the  subroutines  calculating  day  length  and 
instantaneous  irradiance  (ASTRO)  and  instantaneous  gross  assimilation 
(ASSIM),  and  at  three  times  of  the  day  in  the  subroutine  calculating  daily  total 
gross  assimilation  (TOTASS;  Gaussian  integration).  Instantaneous  gross 
assimilation  is  calculated  per  second  and  converted  to  hourly  rates  within 
ASSIM. 

Model  approach  and  organization  are  similar  to  those  used  for  agricultural 
crops  (SUCROSl;  Goudriaan,  Van  Keulen,  and  Van  Laar  1992).  Several 
features  of  a  simulation  model  for  hydrilla  (HYDRIL;  Best  and  Boyd  1996; 

Boyd  and  Best  1996)  and  of  a  general  growth  model  for  submersed  angiosperms 
(SUBANG;  Best  and  Jacobs  1990)  have  been  used. 

MILFO  runs  within  a  FORTRAN  SIMULATION  ENVIRONMENT  (FSE) 
shell.  Version  2.1,  to  enable  easy  handling  of  input  and  output  files  and  rapid 
visualization  of  the  simulation  results  (Van  Kraalingen  1995).  It  can  be  executed 
on  IBM  PC-  ATs  and  compatibles  as  a  stand-alone  version.  Because  of  its 
language  and  simple  structure,  it  will  generally  be  compatible  with  ecosystem 
models  that  accept  FORTRAN. 

The  organization  of  the  model  and  its  subroutines  in  combination  with  the 
FSE  shell  is  illustrated  in  Figure  1. 


Model  Features 

Features  of  MILFO  are  as  follows: 

a.  Phenology  is  tied  indirectly  to  air  temperature  through  development  rate 
and  is,  therefore,  independent  of  day  number;  thus,  the  model  can  be  used 
under  climatological  conditions  ranging  from  temperate  to  tropical. 

b.  Plant  growth  starts  from  the  rhizome/root  crown  system  alone  or  from  the 
same  system  with  wintering  plants. 
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Figure  1 .  Relational  diagram  of  MILFO  and  its  subroutines  in  combination  with 
FSE  shell  (Each  plant  cohort  is  represented  by  a  cohort-specific 
subroutine  (cohort  1  by  MODEL,  cohort  2  by  CHRT2,  and  cohort  3  by 
CHRT3;  only  one  shown),  all  using  same  subroutines  ASTRO, 
TOTASS,  and  ASSIM) 
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c.  Two  plant  cohorts  are  active  in  a  temperate  climate  and  three  cohorts  in 
the  tropics,  depending  on  the  seasonal  input  variables. 

d.  Photosynthetic  response  is  to  instantaneous  irradiance. 

e.  Removal  of  biomass  through  mechanical  harvesting  can  be  calculated. 

/.  Air  or  water  temperatures  must  be  used  to  run  the  model. 

g.  The  model  can  be  used  for  communities  at  various  water  depths,  ranging 
from  0.5  to  6.0  m. 

h.  Plant  parameter  values  and  climatological  variables  can  be  easily 
changed. 
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3  Model  Processes 


Morphology,  Phenological  Cycle, 
and  Development 

Morphology  and  phenological  cycle  of  milfoil 

Eurasian  watermilfoil  is  a  rooted  perennial  with  long,  flexible  stems  and 
finely  dissected  leaves.  The  leaves  are  arranged  in  whorls  around  the  stems. 

The  plant  stems  may  reach  lengths  in  excess  of  4  m  in  summer,  branching  close 
to  the  water  surface  (canopy  formation).  It  has  been  found  in  water  depths  rang¬ 
ing  from  0.2  to  6  m  (Grace  and  Wetzel  1978;  Madsen,  Eichler,  and  Boylen 
1988).  It  occasionally  forms  small  emergent  shoots  from  fragments  starting  on 
the  shore.'  The  current  model  does  not  describe  plants  in  emergent  habit. 

Milfoil  is  able  to  propagate  itself  by  seeds,  by  vegetative  fragmentation,  and 
in  an  evergreen  condition.  Flowering  of  milfoil  in  the  northern  hemisphere 
occurs  from  June  to  November;  one  (Aiken,  Newroth,  and  Wile  1979;  Grainger 
1947;  Carpenter  1980),  two  (Nichols  1971;  Lind  and  Cottam  1969;  Patten  1956), 
and  three  (Grace  and  Wetzel  1978)  flowering  periods  per  year  have  been 
reported.  Flowering  periods  in  southern  areas  have  been  described  as  “less  pre¬ 
dictable”  (Grace  and  Wetzel  1978),  while  they  are  suggested  to  occur  in  the 
tropics  during  the  whole  growth  season  (Zutschi  and  Vass  1973).  Flowering 
usually  coincides  with  peak  biomass  and  is  followed  immediately  by 
autofragmentation/sloughing.  The  production  of  viable  seeds  requires  emersion 
of  the  typically  monoecious  flowering  spikes  (Patten  1954)  with  transfer  of 
pollen  by  wind  as  the  dominant  pollination  mechanism  (Hutchinson  1975). 

Seeds  are  important  in  long-distance  dispersal  and  as  insurance  against  local 
extinction,  but  seed  germination  may  be  delayed  (Guppy  1897;  Patten  1955)  or 
decreased  by  desiccation  (Standifer  and  Madsen  199'^;  seedling  establishment 
appears  to  be  a  particularly  fragile  stage  in  the  life  cycle  (Patten  1956;  Hartleb, 
Madsen,  and  Boylen  1993).  Shoot  fragmentation  is  usually  the  result  of  abscis¬ 
sion  just  after  flowering,  but  it  can  also  be  accidental  (by  boat  contact  or  wave 
action).  Although  shoot  fragmentation  can  be  substantial,  the  number  of 


’  Personal  Communication,  1998,  J.  E.  Titus,  University  of  Binghamton,  New  York. 
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established,  new  plants  originating  from  shoot  fragments  is  relatively  low 
(Madsen  and  Smith  1997).  Fragmentation  is  probably  the  most  important  means 
of  dispersal  within  a  water  body  or  between  nearby  water  bodies.  Milfoil  most 
frequently  winters  in  an  evergreen  form  as  root  crowns  and/or  lower  shoots 
attached  to  the  rhizome  system  (Grace  and  Wetzel  1978;  Madsen,  Eichler,  and 
Boylen  1988;  Madsen  1997)  and  may  maintain  considerable  winter  biomass 
(Stanley  et  al.  1976).  This  species  does  not  form  turions  described  as  important 
hibemacula  of  other  Myriophyllum  species  (M.  exalbescens,  M.  verticillatum, 

M.  heterophyllum;  Grace  and  Wetzel  1978). 


Description  of  deveiopment  and  phenological  cycie  in  MiLFO 

The  phenology  of  a  plant  community,  for  which  development  phase  can  be 
used  as  a  measure,  quantifies  physiological  age  and  is  related  to  its  morphologi¬ 
cal  appearance.  Development  phase  caimot  be  expressed  simply  as  chronologi¬ 
cal  age  because  several  environmental  factors  such  as  temperature  and  stress 
(e.g.,  nutrients,  grazing)  can  speed  up  or  reduce  the  rate  of  phenological  develop¬ 
ment.  Contrary  to  what  is  suggested  by  intuition,  the  rate  of  plant  growth  per  se 
has  no  effect  on  phenological  development,  as  long  as  the  growth  rate  is  not  very 
low  (Penning  de  Vries  et  al.  1989b,  and  citations  therein).  The  concept  of 
deveiopment  phase  is  used  to  characterize  the  whole  plant  community;  it  is  not 
appropriate  for  individual  organs. 

The  response  of  developmental  rate  to  temperature  in  the  current  model  is  in 
accordance  with  the  degree-day  hypothesis  (Thomley  and  Johnson  1990a).  The 
idea  is  as  follows.  The  mean  temperature  r,for  each  day  i  is  measured,  and  a 
sum  h  is  formed  according  to 

h=  s  (7;-r,) 

1=1 

which  includes  only  those  terms  where  is  above  some  threshold  value  . 
When  h  reaches  a  particular  value,  this  signifies  that  a  phase  in  development  is 
complete,  and  this  is  generally  associated  with  a  biological  event  that  occurs 
over  a  short  period  of  time  and  is  readily  observed.  The  day-degree  sum  h 
essentially  integrates  some  underlying  temperature-dependent  processes.  For 
milfoil,  for  example,  there  are  various  phases  in  the  development  of  the  plant, 
and  the  temperature  sum  is  found  to  have  a  certain  value  for  the  successful 
completion  of  each.  The  temperature  threshold  may  be  different  for  each  of 
these  phases.  The  approach  is  based  on  the  notion  of  a  developmental  rate, 
whose  response  to  temperature  is  approximately  linear  over  a  restricted  tempera¬ 
ture  range.  Comparison  with  actual  temperature  responses  found  in  agricultural 
crops  suggests  that  this  is  not  unreasonable,  and  the  method  works  well  in  prac¬ 
tice.  It  is  implicitly  assumed  that  the  organ  possesses  a  developmental  clock  that 
is  proceeding  at  the  rate  kj .  In  general,  it  is  to  be  expected  that  the  development 
rate  may  depend  on  a  number  of  quantities.  This  can  be  represented  by 
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k,  =  f(y,p,E) 


in  which /represents  some  function  of  the  state  variables  V,  parameters  P,  and 
environmental  quantities  E.  The  temperature-sum  rule  works  because  the  most 
important  environmental  variable  is  temperature,  and  the  response  to  temperature 
is  approximately  linear. 

The  phenological  cycle  is  described  using  milfoil  in  Lake  Wingra,  Wisconsin, 
in  1970  as  an  example  (Adams  and  McCracken  1974).  Plant  data  of  this  year 
were  chosen  after  verifying  that  climatological  conditions  did  not  deviate  from 
the  usual  at  that  site. 

Development  phase  (DVS)  is  a  state  variable  in  MILFO.  The  development 
phase  is  dimensionless,  and  its  value  increases  gradually  within  a  growing  sea¬ 
son.  The  development  rate  has  the  dimension  d  \  The  multiple  of  rate  and  time 
period  yields  an  increment  in  phase.  In  the  model,  the  temperature  that  affects 
development  of  milfoil  can  be  chosen  as  equal  to  the  daily  average  air  tempera¬ 
ture  at  the  height  of  the  growing  point  of  the  shoots,  with  a  user-defined  lag- 
period  to  correct  for  deviations  in  temperature  of  the  water  body  in  which  the 
aquatic  community  grows  compared  with  air  temperatures  (7  days  is  nominal).  It 
is  more  accurate  to  use  water  temperatures  for  this  purpose;  but  since  water  tem¬ 
peratures  are  not  always  available  for  the  site  for  which  the  user  wants  to  run  the 
model,  MILFO  can  be  run  using  either  one. 

The  rate  of  phenological  development  can  be  affected  by  temperature  differ¬ 
ently  in  the  vegetative  phase  and  in  the  reproductive  phase.  These  differences 
indicate  that  the  physiological  process  of  development  may  not  be  the  same 
before  and  after  anthesis.  Descriptions  in  literature  of  number  of  flowering  per¬ 
iods  per  year  and  their  timing  in  milfoil  indicate  that  from  June  to  November 
usually  two  flowering  periods,  in  June  and  July,  occur  in  temperate  climates, 
sometimes  three  in  southern  regions,  and  usually  three  in  tropical  climates 
(Zutschi  and  Vass  1973). 

The  following  development  rates  were  derived  from  the  Lake  Wingra  field 
data,  pertaining  to  two  plant  cohorts  each  with  its  own  flowering  period  (Adams 
and  McCracken  1974):  of  0.022  d'^  prior  to  the  first  flowering  period  and  of 
0.015  d'^  subsequently,  at  a  reference  temperature  of  30  °C  and  a  temperature 
threshold  of  3  °C.  These  development  rates  are  considered  as  typical  for  tem¬ 
perate  regions. 

For  milfoil  populations  in  the  tropics,  the  same  development  rates  and  tim¬ 
ings  as  in  temperate  regions  were  applicable,  but  a  third  plant  cohort  had  to  be 
added  to  accommodate  the  third  flowering  period  and  usually  high  August  bio¬ 
mass  in  India  (Zutschi  and  Vass  1973).  The  milfoil  development  rates  are  some¬ 
what  higher  than  those  found  for  hydrilla  (0.012  d'^  at  the  same  reference 
temperature  and  threshold  temperature  as  used  for  milfoil). 

The  development  phase  has  the  value  zero  when  the  simulation  starts  at  the 
first  Julian  day  number  (Tables  1  and  2).  The  simulation  starts  using  observed 
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Table  1 

Relationship  Between  Development  Phase  (DVS)  of  Milfoil,  Day  of  Year,  and  3  °C  Day- 
Degree  Sum  for  a  Temperate  Climate  (DVRVT=  0.022;  DVRRT=  0.015) 

Developmental  phase 

Day  Number 

3  ®C  Day-Degree 
Sum 

Description 

DVS  Value 

First  Julian  day  number  ->  sprouting,  initiation  elongation,  and 
leaf  expansion  COHORT1 

0  ->  0.375 

0->114 

1  ->191 

Sprouting,  initiation  elongation,  and  leaf  expansion  ->  floral 
initiation,  anthesis,  and  induction  of  senescence  COHORT1 

0.376 ->1.000 

115->16S 

19S  ->  900 

Floral  initiation,  anthesis,  and  induction  of  senescence 
->  senescence  COHORT1 

1.001  ->1.630 

163  ->  21S 

901  ->  S01S 

Senescence  ->  senesced  COHORT  1 

1.631->S.000 

SI  3  ->  S45 

S013->  2669 

Sprouting,  initiation  elongation,  and  leaf  expansion  ->  floral 
initiation,  anthesis,  and  induction  of  senescence  COHORT2 

1.001 ->1.630 

163  ->  SIS 

901 ->  2012 

Floral  initiation,  anthesis,  and  induction  of  senescence 
->  anthesis  and  senescence  COHORTS 

1.631  ->  S.OOO 

S13->  S45 

2013 ->  2669 

Senescence  ->  senesced  COHORTS 

S.001  ->  S.570 

S46  ->  365 

2670  ->  3508 

Senesced  COHORT  1  and  S 

S.570 

365 

3508 

Note:  Calibration  was  on  field  data  on  biomass  ( Adams  and  McCracken  1974)  and  on  water  transparency,  temperature,  and 
irradiance  from  Lake  Wingra,  Wl,  1970  (Lee  eind  Kluesener  1972). 

weights  of  plants  and  rhizome/root  crowns  as  initial  values.  Initial  plant  weights 
have  been  set  equal  to  the  observed  shoot  weight  early  in  spring,  which  is 
believed  to  give  a  fair  approximation.  Since  the  initial  weight  of  the  rhizome/ 
root  crown  system  had  not  been  measured  in  the  calibration  data  set,  this  weight 
has  been  set  equal  to  50  g  DW  m'^  found  for  a  similar  milfoil  community  in  the 
same  lake  in  1977  (Smith  and  Adams  1986).  The  rhizome/root  crown  system  is 
the  common  basis  from  which  milfoil  plant  cohorts  develop.  Plant  cohorts  are 
plant  groups  exhibiting  the  same  phenological  cycle,  and  plants  are  considered  as 
units  composed  of  roots,  stems,  and  leaves,  excluding  the  rhizome/root  crown 
system.  If  simulation  of  the  community  at  another  site  is  desired,  the  simulation 
can  start  from  other  initial  biomass  values,  either  from  the  rhizome/root  crown 
system  only  or  with  wintering  plant  biomass  present. 

For  a  milfoil  community  in  a  temperate  climate  (Table  1),  the  sprouting  of  the 
rhizome/root  crown  system,  i.e.,  the  initiation  of  growth  activity,  occurs  at  DVS 
0.375.  Sprouts  of  plant  cohort  1  develop  through  remobilization  of  carbohy¬ 
drates  from  the  rhizome/root  crown  system.  The  sprouts  elongate  rapidly  to  the 
water  surface  and  form  a  canopy  in  the  upper-water  layers.  Anthesis  of  cohort  1 
is  initiated  at  DVS  1.000  and  finishes  at  DVS  1.630,  just  before  downward  car¬ 
bohydrate  translocation  and  senescence  are  initiated.  Translocation  and  senes¬ 
cence  of  cohort  1  set  in  at  DVS  1.631  and  continue  until  DVS  2.000.  Sprouting 
of  cohort  2  starts  when  translocation  and  scenescence  of  cohort  1  have  set  in. 
This  timing  is  based  on  the  assumption  that  at  that  time,  apical  dominance  by  the 
existing,  senescing  shoots  is  broken  and,  consequently,  new  shoots  can  develop. 
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Table  2 

Relationship  Between  Development  Phase  (DVS)  of  Milfoil,  C 
Degree  Sum  for  a  Tropical  Climate  (DVRVT=  0.022;  DVRRT= 

lay  of  Year,  an 
0.015) 

Id  3  °C  Day- 

Developmental  phase 

Description 

DVS  Value 

Day  Number 

3  "C  Day-Degree 
Sum 

First  Julian  day  number  ->  sprouting,  initiation  elongation,  and 
leaf  expansion  COHORT1 

0  ->  0.S75 

0->  25 

1  ->  4S1 

Sprouting,  initiation  elongation,  and  leaf  expansion  ->  floral 
Initiation,  anthesis,  and  induction  of  senescence  COHORT1 

0.S76  ->  1.000 

26  ->  61 

4S2  ->  116S 

Floral  initiation,  anthesis,  and  induction  of  senescence 
->  senescence  COHORT1 

1.001  ->1.6S0 

62  ->  162 

1164 ->  S844 

Senescence  ->  senesced  COHORT  1 

1.6S1  ->  2.000 

16S  ->188 

3845  ->  4490 

Sprouting,  initiation  elongation,  and  leaf  expansion  ->  floral 
initiation,  anthesis,  and  induction  of  senescence  COHORT2 

1.001  ->1.6S0 

62 ->162 

1164->  S844 

Floral  initiation,  anthesis,  and  induction  of  senescence 
->  anthesis  and  senescence  COHORT2 

1.6S1  ->2.000 

16S  ->188 

S845  ->  4490 

Senescence  ->  senesced  COHORT2 

2.001  ->  2.570 

164  ->  2SS 

4491  ->  5492 

Sprouting,  initiation  elongation,  and  leaf  expansion  ->  floral 
initiation,  anthesis,  and  Induction  of  senescence  COHORTS 

2.001  ->  2.447 

164  ->  22S 

4491  ->  527S 

Floral  initiation,  anthesis,  and  induction  of  senescence 
->  senescence  COHORTS 

2.448  ->  S.500 

224  ->  S07 

5274 ->7125 

Senescence  ->  senesced  COHORTS 

S.501  ->  4.141 

S08  ->  S65 

7126  ->  8254 

Senesced  COHORT  1 ,2,  and  S 

4.141 

S65 

8254 

Note:  Calibration  was  on  field  data  on  biomass  from  Kashmir  lakes,  India,  1970s  (Zutschi  and  Vass  197S)  and  climatological  data 
from  Patancheru,  India,  1978. 

Sprouting  of  cohort  2  occurs  from  growing  points  on  the  rhizome/root  crown 
system.  Anthesis  of  cohort  2  is  initiated  at  DVS  1.631  and  finishes  at  DVS 
2.000.  Translocation  and  senescence  of  cohort  2  set  in  at  DVS  2.001  and  con¬ 
tinue  until  the  end  of  the  year. 

For  a  milfoil  community  in  the  tropics  (Table  2),  it  proved  impossible  to  gen¬ 
erate  the  high  levels  of  shoot  and  rhizome/root  crown  biomass  reported  (Zutschi 
and  Vass  1973)  with  two  plant  cohorts  active  since  the  second  plant  cohort  had 
already  senesced  in  May.  However,  proper  biomass  levels  and  timing  were 
attained  with  three  plant  cohorts  active,  the  third  cohort  being  switched  on  at 
latitudes  less  than  33  °N.  It  is  possible  that  a  particular  plant  process,  like 
sprouting,  is  sensitive  to  day  length  and  that  this  process  decides  for  the  popula¬ 
tion  to  activate  another  cohort.  However,  since  the  authors  are  not  aware  of  pub¬ 
lications  on  this  topic  for  milfoil,  the  switch  has  been  set  at  the  cut-off  latitude 
for  tropical  areas.  Plant  cohorts  in  tropical  regions  behave  similar  in  terms  of 
DVS  to  those  in  temperate  regions,  except  that  tropical  cohorts  require  on  aver¬ 
age  a  1.6  X  higher  3°  degree-day  sum  to  complete  their  individual  life  cycle  than 
temperate  cohorts. 
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Maximum  Biomass  and  Piant  Density 


Seasonal  biomass  maxima  have  been  reported  to  vary  considerably  over  time 
and  space.  In  temperate  climates,  sometimes  one,  but  usually  two,  biomass 
peak(s)  were  found  per  growth  season.  Biomass  maxima  appear  to  be  related  to 
flowering  period.  One  distinct  biomass  maximum  has  been  reported  for  tropical 
areas  (India),  while  flowering  started  in  May  and  continued  during  the  growth 
season.  The  highest  standing  crop  of  2,283  g  DW  m'^  has  been  found  in  Fish 
Lake,  Wisconsin  (Budd,  Lillie,  and  Rasmussen  1995),  and  similar  values  have 
been  reported  for  the  more  southern  Lake  Guntersville,  Alabama,  in  1972 
(Stanley  et  al.  1976).  This  maximum  biomass  value  found  has  been  used  to  form 
the  upper  limit  of  plant  biomass  in  the  model. 

Generally,  biomass  production  of  milfoil  is  far  more  constrained  by  plant- 
inherent  factors,  light  and  space  availability  and  temperature,  than  by  plant  dens¬ 
ity.  As  the  season  progresses,  the  individual  plant  size  increases  along  with  the 
areal  biomass,  and  thinning  of  shoots  caused  by  intraspecific  interference  results 
in  an  inverse  relationship  between  plant  size  and  plant  density  (Lind  and  Cottam 
1969). 

However,  since  initial  plant  density  is  required  as  an  input  variable  into  the 
model,  a  feasible  plant  density  under  field  conditions  had  to  be  found.  A  range 
of  3  to  32  so-called  “plant  clumps”  m  ^  consisting  of  a  variable  number  of  stems, 
were  determined  for  a  milfoil  community  in  Fish  Lake,  Wisconsin,  in  the  sum¬ 
mers  of  1990-92  (Budd,  Lillie,  and  Rasmussen  1995).  The  mean  value  of 
11  plants  m with  clump  used  synonymously  to  plant,  has  been  used  in  the 
model. 

In  MILFO,  plant  density  has  been  set  to  11  plants  m'^.  This  implies  that  plant 
density  is  always  11  m'^  at  the  beginning  of  the  growth  season,  and  that  biomass 
is  redistributed  over  11  plants  m'^  if  wintering  plants  are  present. 

Wintering  and  Sprouting  of  Rhizomes/ 

Root  Crowns  and  Growth  of  Sprouts 
to  Water  Surface 

Rhizome/root  crown  tissues  were  the  main  storage  area  for  carbohydrates  in 
wintering  milfoil.  Starch  concentrations  may  reach  20  percent,  with  total  non- 
structural  carbohydrates  (TNC)  concentrations  of  up  to  30-40  percent  (Titus  and 
Adams  1979b;  Madsen  1997).  Rhizome/root  crown  biomass  tended  to  be  higher 
in  spring  and  in  autumn  than  during  the  rest  of  the  year  and  showed  an  inverse 
relationship  with  plant  cohort  biomass.  It  fluctuated  between  12  and  400  g  DW 
m'^  because  of  seasonal  changes  (Madsen  1997).  Rhizome/root  crown  biomass  of 
the  milfoil  community  in  Lake  Wingra  amounted  to  50  g  m‘^  in  winter  1977 
(Smith  and  Adams  1986),  while  it  was  relatively  constant  in  tropical  regions, 
varying  between  32  and  48  g  DW  m‘^  (Zutschi  and  Vass  1973), 
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In  the  model,  rhizome/root  crown  weight  decreases  by  sprouting  of  growing 
points,  which  transform  into  plants,  by  respiration,  by  a  plant-inherent  sloughing 
process,  and,  possibly,  by  grazing  by  waterfowl  or  other  organisms,  and  it 
increases  by  downward  carbohydrate  translocation. 

Sprouting  or  regrowth  potential  of  the  rhizome/root  crown  system  is  usually 
high  and  occurs  early  in  the  season.  Sprouting  in  southern  areas  like  Texas  (lati¬ 
tude  33  °N,  longitude  97  °E)  has  been  reported  to  occur  already  in  March  (Mad¬ 
sen  1997).  In  northern  areas,  the  timing  of  sprouting  may  be  similar,  but  no 
observations  confirming  this  have  been  made  (or  published)  probably  since  at 
that  time,  water  temperatures  are  still  very  low,  impeding  field  work.  Actual 
sprouting  frequency  under  natural  conditions  is  unknown.  Sprouting  frequency 
in  an  established  community  is  probably  not  important,  as  long  as  the  final  plant 
density  of  11  plants  m'^  is  somehow  reached,  since  plant  density  tends  to  play  a 
lesser  role  in  biomass  production  compared  with  space  availability  (see  Maxi¬ 
mum  Biomass  and  Plant  Density).  Sensitivity  to  day  length  at  which  the 
rhizome/root  crown  systems  sprout,  or  triggering  by  red-far  red  ratio,  has  not 
been  reported. 

It  is  to  be  expected  that  the  rhizome/root  crown  system  requires  continuous 
maintenance,  but  that  maintenance  processes  proceed  at  a  low  level  of  activity 
because  of  the  relatively  high  carbohydrate  concentrations  that  are  cheap  in 
maintenance  costs  (Penning  de  Vries  and  Van  Laar  1982b). 

Sloughing  or  death  rates  of  rhizomes/root  crowns  have  not  been  published  so 
far.  A  death  rate  value  has  been  derived  from  observations  on  terrestrial  rhizome 
systems  where  annual  turnover  rates  were  found  to  be  approximately  four  times 
less  than  those  of  aboveground  plant  biomass  in  the  growth  season,  but  could 
drop  with  a  factor  of  1/100  in  inactive  periods  (Vogt,  Vogt,  and  Bloomfield 
1991).  Following  this  approach,  a  tentative  relative  death  rate  of  0.00042  d'‘  was 
calculated  (g  DW  g  DW  ‘  d  ‘),  being  1/100  of  the  plant  death  rate.  The  latter 
value  is  far  lower  than  that  of  0.36  d'^  estimated  for  hydrilla  tubers  from  simula¬ 
tions  (Best  and  Boyd  1996).  However,  the  death  rate  of  hydrilla  tubers  may  be 
an  overestimate  since  death  by  grazing  and/or  disturbance  of  sediments  was 
included  in  that  overall  death  rate,  and  grazing  of  tubers,  e.g.,  by  waterfowl,  is 
usually  high  (Jupp  and  Spence  1977;  Scheffer,  Bakema,  and  Wortelboer  1993). 
Effects  of  grazing  on  the  milfoil  rhizome/root  crown  systems  are  unknown,  but 
expected  to  be  far  lower  than  on  hydrilla  tubers. 

In  MILFO,  initial  rhizome/root  crown  biomass  has  been  set  at  50  g 
AFDW  m'^,  equal  to  the  below-ground  biomass  measured  at  1.5-m  rooting  depth 
in  Lake  Wingra  in  1977  (Smith  and  Adams  1986)  and  equal  to  the  lowest  shoot 
biomass  found  in  1970  (Adams  and  McCracken  1974).  Sprouting  is  a  function 
of  devel-opment  phase  through  the  3  °C  day-degree  sum;  it  occurs  between  DVS 
0.375  and  the  flowering  period  for  cohort  1,  between  DVS  1.001  and  the  flower¬ 
ing  period  for  cohort  2,  and,  when  active,  between  DVS  2.001  and  the  flowering 
period  for  cohort  3.  Sprouting  frequency  has  been  set  equal  to  the  number  of 
plants  per  surface  area,  i.e.,  at  11  sprouts  m'^  (sprout  is  used  here  synonymously 
with  plant  clump). 
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Remobilization  can  occur  provided  the  weight  of  the  rhizome/root  crown 
system  is  greater  than  the  critical  rhizome  weight.  The  critical  rhizome  weight 
value  is  the  lowest  value  published,  i.e.,  12.0  g  DW  m'^  (Madsen  1997). 

The  rhizomes  sprout  by  remobilization,  i.e.,  conversion  of  part  of  their  carbo¬ 
hydrate  reserves  into  sprout  material,  via  a  relative  conversion  rate  of  rhizome- 
to-plant  (ROC),  with  the  same  value  as  derived  for  conversion  of  hydrilla  tubers 
(0.0576  g  CH2O  g  rhizome/root  crown-DW'^  d  ' ;  for  calculation,  see  Best  and 
Boyd  1996).  These  carbohydrates  are  allocated  to  the  plant  organs  according  to 
a  fixed  biomass  allocation  pattern  (see  next  section).  It  is  assumed  that  the 
sprouts  can  elongate  up  to  the  water  surface  by  mere  remobilization  processes, 
not  even  requiring  photosynthetic  products,  since  potential  sprout  elongation  has 
been  estimated  to  be  12  m  sprout-length  g  DW"^  (for  hydrilla,  cf.  Best  and  Boyd 
1996;  it  is  assumed  to  be  similar  for  milfoil). 

After  reaching  the  water  surface,  canopy  formation  takes  place  and  photosyn¬ 
thesis  proceeds. 

Maintenance  processes  are  treated  in  the  next  section. 

The  relative  rhizome  death  rate  has  been  set  at  0.00042  d'^  (on  dry  weight 
basis). 

A  relational  diagram  illustrating  the  wintering  and  sprouting  rhizomes/root 
crowns  of  milfoil  is  shown  in  Figure  2. 

-TRANSl  +  REMOBl  +  MAINRT 

TWGRIZ  =  INTGRL(TWGRIZ,  -  [  (- . -)  +  (RDRIZ  x  TGRIZ)],  BELT) 

1.242 

IF  (DVS.  GE.  0.376.  AND.  DVS.  LT.  1.0)  THEN 

IF  (TWGRIZ.  GT.  CRRIZ)  THEN 

REMOBl  =  ROC  X  TWGRIZ 

TWRIZD  =  INTGRL  (TWRIZD,  RDRIZ,  BELT) 

where 

TWGRIZ  =  total  live  dry  weight  of  rhizome/root  crown  system  of  current  day 
(gDWm-2) 

REMOBl  =  remobilization  rate  of  carbohydrates  cohort  1  (g  CH2O  m'^  d'^) 

MAINRT  =  maintenance  respiration  rate  of  rhizome/root  crown  system 
(gCH20m-M-i) 

1.424  =  assimilate  requirement  for  rhizome  dry  matter  production 
(g  CH2O  g  DW‘^;  see  Appendix  C) 

RDRIZ  =  relative  death  rate  of  rhizome/root  crown  system  (d’^) 
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Figure  2.  Relational  diagram  illustrating  wintering  and  sprouting  of 
rhizomes/root  crowns  in  milfoil 


Chapter  3  Model  Processes 


TGRIZ  =  total  live  dry  weight  rhizome/root  crown  system  of  previous  day 
(g  DW 

DVS  =  developmental  phase  of  plant  (-) 

CRRIZ  =  critical  weight  of  rhizome/root  crown  system  (g  DW  m'^ 

ROC  =  relative  conversion  rate  of  rhizome/root  crown  into  plant  material 
(gCH,OgDW-‘d‘) 

TWRIZD  =  total  weight  of  dead  rhizome/root  crown  system  (g  DW  m'^ 

Light,  Photosynthesis,  Maintenance,  Growth, 
and  Assimiiate  Partitioning  in  Miifoii  Piants 

Light 

The  measured  daily  total  irradiance  (wavelength  300-3,000  nm)  is  used  as 
input  for  the  model.  Only  half  of  the  irradiance  reaching  the  water  surface  is 
photosynthetically  active  and  is  therefore  used  to  calculate  CO  2  assimilation.  Six 
percent  of  the  irradiance  is  reflected  by  the  water  surface  (Golterman  1975). 

The  subsurface  irradiance  is  attenuated  by  dissolved  substances  and  particles 
within  the  water  column  resulting  in  a  site-  and  season-specific  extinction  coeffi¬ 
cient.  Moreover,  the  vertical  profiles  of  the  radiation  within  the  community  lay¬ 
ers  are  characterized.  The  absorbed  irradiance  for  each  horizontal  community 
layer  is  derived  from  these  profiles.  The  community-specific  extinction  coeffi¬ 
cient  K  is  assumed  to  be  constant  throughout  the  year  and  given  a  value  of 
0.006  m^  g  DW  *  measured  in  the  milfoil  community  in  Lake  Wingra  (Titus  and 
Adams  1979a).  Another  higher,  community-specific  extinction  coefficient  of 
0.01  m^  g  DW  *  has  been  published  by  Ikusima  (1970)  for  a  milfoil  community  in 
Japan,  which  may  indicate  that  plants  at  lower  latitudes  have  thinner  leaves. 

The  incoming  irradiance  is  attenuated  by  the  shoots,  part  of  which  is 
absorbed  by  the  photosynthetic  plant  organs,  i.e.,  the  leaves. 
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where 


IRZ(i)  =  photosynthetic  active  part  of  total  irradiance  on  top  of  depth  layer 
i  (J  s'^) 

TL  =  thickness  depth  layer  (0.10  m) 

L  =  light  extinction  coefficient  of  water  (m'‘) 

K  =  plant-specific  extinction  coefficient  (m^  g  DW  *^) 

SC  =  shoot  matter  (g  DW  per  0.1  m  stratum  of  a  m^  water  column) 
lABS(i)  =  total  irradiance  absorbed  in  depth  layer  i  (J  m’^  s  ‘) 
lABSL(i)  =  total  irradiance  absorbed  by  plant  shoots  in  depth  layer  i  (J  m'^s  ’) 
FL  =  leaf  dry  matter  allocation  to  each  layer  of  the  plant  (relative) 


Photosynthesis 

The  instantaneous  rates  of  gross  assimilation  are  calculated  from  the 
absorbed  light  energy  and  the  photosynthesis  light  response  of  individual  shoot 
apices,  here  used  synonymously  to  leaves. 

The  photosynthesis  light  response  of  leaves  is  described  by  the  exponential 
function 

-EE  X  lABS  X  3600 

FGL  =  SC,  X  AMAX  x  (1  -  exp  [ - - ]) 


where 


FGL  =  gross  assimilation  rate  per  depth  layer  (g  CO2  m'^  h'^) 

SC(i)  =  standing  crop  in  depth  layer  i  (g  DW  m'^  layer  ■^) 

AMAX  =  actual  CO2  assimilation  rate  at  light  saturation  for  individual  shoots 
(gC02gDW-ih-‘) 

EE  =  initial  light-use  efficiency  for  shoots  (g  CO2  absorbed) 

For  photosynthetic  activity  at  light  saturation  (AMAX),  the  value  of  0.0165  g 
CO2  g  DW*  h  *  was  used.  This  value  is  equal  to  the  field  AMAX  measured  in 
Lake  Wingra  in  May  1971,  at  pH  8  and  a  total  alkalinity  of  190  mg  L  *  (Adams 
and  McCracken  1974).  It  is  slightly  higher  than  field  values  measured  for 
hydrilla  in  water  in  equilibrium  with  atmospheric  CO2  (0.0158  g  CO2  g  DW  *  h  *; 
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Bowes,  Holaday,  and  Haller  1979;  Van,  Haller,  and  Bowes  1976).  Light-  and 
carbon-saturated  photosynthetic  rates  can  be  far  higher  (Van,  Haller,  and  Bowes 
1976),  suggesting  that  photosynthetic  activity  in  lakes  like  Lake  Wingra,  where 
Die  concentrations  are  in  the  range  of  0.8  to  3.5  mmol  with  a  pH  of  7.6  to  9.4 
(Lee  and  Kluesener  1972),  can  be  carbon  limited. 

For  photosynthetic  light-use  efficiency  (EE),  a  value  of  11.10  *  g  COj  J  *, 
typical  for  C3  plants,  was  used  (Penning  de  Vries  and  Van  Laar  1982a).  Substi¬ 
tuting  the  appropriate  value  for  the  absorbed  photosynthetically  active  radiation 
yields  the  assimilation  rate  for  each  specific  shoot  layer. 

Gross  assimilation  rate  at  light  saturation  shows  a  distinct  seasonal  pattern 
and  tends  to  decrease  with  aging  (Adams  and  McCracken  1974;  Adams,  Titus, 
and  McCracken  1974).  Although  a  function  describing  this  relationship 
(AMDVST)  has  been  included  in  the  model,  it  is  not  active  in  the  nominal  ver¬ 
sion  (it  has  the  value  of  1)  since  by  running  the  model  it  turned  out  not  to  be 
quantitatively  important.  Gross  assimilation  in  milfoil  tends  to  decrease  from 
apex  to  stem  base  (Adams,  Titus,  and  McCracken  1974).  A  function  describing 
this  relationship  (^DFT)  has  been  included  in  the  model,  but  is  not  active  in 
the  nominal  version  (it  has  the  value  of  1)  since  it  also  turned  out  not  to  be  quan¬ 
titatively  important. 

A  reduction  factor,  REDAM,  can  be  used  to  take  the  effects  of  daily  changes 
in  pH  and  oxygen  concentrations  on  AMX  into  account,  by  reducing  the  AMX 
by  a  factor  between  0  and  1  for  the  whole  day.  REDAM  currently  has  the  value 
of  0.5  since  it  appears  that  pH  in  milfoil  communities  in  Lake  Wingra  usually 
oscillates  around  8.8  (Adams  and  McCracken  1974),  causing  a  50-percent  reduc¬ 
tion  in  photosynthetic  activity  (Titus  and  Stone  1982).  Milfoil  appears  to  be 
relatively  insensitive  to  changes  in  oxygen  concentration  (a  reduction  in  the  net 
photosynthetic  rate  of  only  5  percent  was  observed  because  of  a  change  in  oxy¬ 
gen  concentration  from  1  to  21  percent  at  15  //m  CO2;  Van,  Haller,  and  Bowes 
1976). 

A  fitted,  relative  function,  AMTMPT,  describes  the  effect  of  daytime  tem¬ 
perature  on  photosynthetic  activity.  This  function  has  its  optimum  at  35  °C  and 
is  based  on  the  photosynthetic  response  of  milfoil  to  temperature  (Titus  and 
Adams  1979a;  confirmed  by  Stanley  and  Nailor  1972;  see  Appendix  C). 

The  instantaneous  rate  of  gross  assimilation  over  the  height  of  the  community 
is  calculated  by  relating  the  assimilation  rate  per  layer  to  the  community-specific 
biomass  distribution  and  by  subsequent  integration  of  all  individually  0.1-m-high 
community  layers. 

The  daily  rate  of  gross  assimilation  is  calculated  by  using  the  Gaussian  inte¬ 
gration  method.  This  method  specifies  the  discrete  points  at  which  the  value  of 
the  function  to  be  integrated  has  to  be  calculated  and  the  weighting  factors  that 
must  be  applied  to  these  values  to  attain  minimum  deviation  from  the  analytical 
solution.  A  three-point  method  performs  very  well  for  calculating  daily  total 
assimilation  (Goudriaan  1986;  Spitters  1986). 
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Maintenance,  growth,  and  assimilate  partitioning 

Maintenance.  Some  of  the  carbohydrates  formed  are  respired  to  provide 
energy  for  maintaining  the  existing  plant  components.  The  maintenance  costs 
increase  with  metabolic  activity,  probably  because  of  higher  enzyme  turnover 
and  higher  transport  costs  (Penning  de  Vries  1975). 

The  maintenance  cost  can  be  estimated  from  the  chemical  composition  of  the 
plant.  Typical  maintenance  coefficients  for  various  plant  organs  have  been 
derived,  based  on  numerous  chemical  determinations  in  agricultural  crops.  They 
typically  range  from  0.010  to  0.016  g  CH2O  g  AFDW^  d'^  (Penning  de  Vries  and 
Van  Laar  1982b). 

In  MILFO,  the  maintenance  coefficients  mentioned  above  are  used  to  calcu¬ 
late  the  maintenance  requirement  of  the  plant  cohorts.  A  lower  maintenance 
coefficient  of  0.005  g  CH2O  g  AFDW^  is  used  for  the  rhizome/root  crown  sys¬ 
tem,  considered  to  be  similar  in  respiration  to  stems  with  coefficients  <0.007 
(Penning  de  Vries  et  al.  1989a). 

Higher  temperatures  expedite  turnover  rates  of  plant  tissues  and  increase 
maintenance  costs.  A  temperature  increase  of  10  °C  usually  increases 
maintenance  respiration  by  a  factor  of  about  2  up  to  temperatures  that  usually 
kill  plants  (45  to  60  °C;  Qio=2  at  a  reference  temperature  30  °C;  Penning  de 
Vries  et  al.  1989a). 

Maintenance  respiration  in  MILFO  has  been  related  to  temperature  by  a  fac¬ 
tor,  TEFF,  which  has  the  value  of  1  between  5  and  20  °C  (increases  twofold  with 
every  10  °C  above  a  reference  temperature  of  20  "C  (Thomley  and  Johnson 
1990a)  and  increases  linearly  from  0.0001  to  1  between  0  and  5  °C).  The  value 
of  2  appears  to  be  a  reasonable  average,  but  lower  and  higher  Qjo  values  have 
been  reported  also  (Amthor  1984).  The  currently  used  Qk,  value  is  lower  in  the  0 
to  20  "C  range  than  2.28  found  for  a  Qio  of  dark  respiration  in  milfoil  (Grace  and 
Wetzel  1978);  however,  the  latter  process  includes  growth  processes. 

Equations  describing  maintenance  costs  for  milfoil  plant  cohorts  (1, 2,  or  3)  are: 

MAINTS  =  0.016  X  TWLG  +  0.010  x  TWSG  +  0.015  x  TWRG 

MAINT  =  MAINTS  x  TEFF 
where 

MAINTS  =  maintenance  respiration  rate  plant  at  reference  temperature 
(gCH20m-2d-‘) 

TWLG  =  total  dry  weight  of  live  leaves  (g  DW  m'^ 
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TWSG  =  total  dry  weight  of  live  stems  (g  DW  m'^ 

TWRG  =  total  dry  weight  of  roots  (g  DW 
Equations  describing  maintenance  costs  for  the  rhizome/root  crown  system  are: 

MAINRT  =  0.005  x  TWGRIZ  x  TEFF 
where 

MAINRT  =  maintenance  respiration  rate  rhizome/root  crown  system  at 
reference  temperature  (g  CH2O  m'^  d  ‘) 

TWGRIZ  =  total  dry  weight  of  rhizome/root  crown  system  of  current  day 
(gDWm-^) 

Growth.  Assimilates  in  excess  of  maintenance  costs  are  available  for  con¬ 
version  into  structural  plant  material.  In  this  conversion  process  of  the  glucose 
molecule,  COj  and  HjO  are  released.  The  assimilates  required  to  produce  one 
unit  weight  of  any  particular  plant  organ  can  be  calculated  from  its  chemical 
composition  and  the  assimilate  requirements  of  the  various  chemical  compo¬ 
nents.  Typical  values  are  1.46  g  CH2O  g  DW'^  for  leaves,  1.51  for  stems,  and 
1.44  for  roots  (Penning  de  Vries  and  Van  Laar  1982b;  Penning  de  Vries  et  al. 
1989a),  confirmed  by  Griffin  (1994).  At  higher  temperatures,  the  conversion 
processes  are  accelerated,  but  the  pathways  are  identical.  The  recently  deter¬ 
mined  construction  costs  for  several  submersed  plant  species,  using  a  different 
method  (Williams  et  al.  1987),  are  generally  lower,  ranging  jfrom  0.99  to  1.11 
(Spencer,  Ryan,  and  Ksander  1997).  However,  the  latter  plants  appear  to  be 
relatively  poor  in  nitrogen,  and  transport  costs  have  not  been  included,  both  fac¬ 
tors  that  may  have  contributed  to  this  lower  cost  calculated. 

In  MILFO,  the  construction  costs  typical  for  agricultural  plants  have  been 
used  since  construction  costs  calculated  for  milfoil  shoots  with  an  average  chem¬ 
ical  composition  were  similar  to  those  in  agricultural  plants,  i.e.,  1.54  CH2O  g 
DW"'  (see  Appendix  C). 

The  following  equation  describes  growth: 

QTW  =  ^  ^  GPHOT  -  TRANSl  -  MAINT) 

ASRQ 

where 


GTW  =  dry  matter  growth  rate  of  vegetation  (plants  excluding 
rhizome/root  crown  system;  g  DW  m"^  d"') 

GPHOT  =  daily  total  gross  assimilation  rate  of  community  (g  CH2O  m"^  d"') 

REMOB  1  =  remobilization  rate  of  carbohydrates  cohort  1  (g  CH2O  m'^  d"') 
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TRANSl  =  translocation  rate  of  carbohydrates  cohort  1  (g  CH2O  d 
MAINT  =  maintenance  respiration  rate  of  vegetation  (g  CH2O  m'^  d'^) 

ASRQ  =  assimilate  requirement  for  plant  dry  matter  production  (g  CHjO 
gDW-i) 

Assimilate  partitioning.  Assimilate  partitioning  is  the  process  by  which 
assimilates  available  for  growth  are  allocated  to  leaves,  stems,  roots,  and/or  stor¬ 
age  organs.  The  distribution  pattern  is  a  function  of  physiological  age. 

In  MILFO,  the  assimilate  allocation  pattern  has  been  used  synonymously  with 
the  biomass  allocation  pattern.  This  pattern  is  assumed  to  be  followed  after  the 
shoot  tips  have  reached  the  water  surface  and  not  to  change  with  physiological 
age  (only  summer  values  on  biomass  partitioning  were  available).  The  assimi¬ 
late  allocation  has  been  set  at  0.47  of  total  net  growth  (excluding  rhizome/root 
crown  system)  to  leaves,  0.47  to  stems,  and  0.06  to  roots.  These  values  have 
been  derived  from  the  compartmentalization  of  biomass  over  plant  organs  in  a 
well-developed  milfoil  community,  with  shoots  composed  of  50  percent  by 
leaves  and  50  percent  by  stems  (Budd,  Lillie,  and  Rasmussen  1995).  A  contribu¬ 
tion  of  0.06  to  total  biomass  by  roots  was  chosen  since  no  information  on  the 
roots  of  the  same  vegetation  was  given,  but  root  biomass  is  known  to  be  usually 
small  (similar  to  the  contribution  of  roots  to  total  plant  biomass  in  hydrilla;  Best 
and  Boyd  1996).  Contributions  of  leaves  and  stems  to  total  biomass  were  recal¬ 
culated  proportionally. 

The  following  equation  describes  biomass  allocation  to  plant  organs: 

GRT  =  FRTxGTW 

GST  =  FSTxGTW 

GLV  =  FLVxGTW 
where 

GRT,  GST,  and  GLV  =  dry  matter  growth  rates  of  roots,  stems,  and  leaves, 
respectively  (g  DW  m'^  d'^) 

FRT,  FST,  and  FLV  =  fraction  of  total  dry  matter  allocated  to  roots,  stems, 
and  leaves,  respectively  (relative) 

GTW  =  dry  matter  growth  rate  of  the  vegetation  (plants 

excluding  rhizome/root  crown  system;  g  DW  m'^  d'^) 

In  adolescent  milfoil  plants,  shoot  biomass  is  usually  present  for  61  percent  in 
the  upper  0.5-m  water  column,  distributed  for  10  percent  in  the  upper  0.1-m 
layer,  for  16  and  17  percent  in  both  successive  layers  2  and  3,  and  for  10  and 
8  percent  in  both  successive  layers  4  and  5  (Adams,  Titus,  and  McCracken 
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1974).  These  values  form  the  basis  for  the  dry  matter  allocation  per  depth  layer 
over  the  vertical  axis,  from  water  surface  to  0.5-m  depth.  The  values  of  this 
function  (DMPC)  are  read  from  the  input  file  and  can  be  changed  by  the  user. 
Dry  matter  allocation  to  the  lower  water  layers  is  equal  up  to  a  total  biomass 
share  of  5  percent.  The  remaining  biomass  is  divided  proportionally  over  all 
water  layers.  Vertical  biomass  distribution  pattern  is  recalculated  and  redistri¬ 
buted  by  MILFO  when  a  rooting  depth  other  than  the  nominal  one  (1.5  m)  is 
chosen. 

A  relational  diagram  illustrating  photosynthesis,  respiration,  and  biomass 
formation  of  milfoil  is  shown  in  Figure  3. 


Flowering,  Translocation,  and  Senescence 

The  occurrence  of  flowering  affects  subsequent  metabolic  activity  of  the 
vegetation.  The  timing  of  flowering  is,  therefore,  extremely  important  for  the 
physiological  activity  and  biomass  formation,  while  the  actual  investment  of  dry 
matter  in  flowers  and  seeds  proves  to  be  only  minor  (Madsen  1997). 

After  flowering,  scenescence  sets  in,  and  a  considerable  part  of  net  produc¬ 
tion  is  translocated  downwards  to  the  rhizome/root  crown  system,  while  the 
remainder  of  net  production  is  allocated  according  to  the  above-mentioned  key. 

The  translocated  material  consists  mainly  of  carbohydrates  and  proteins  and 
is  largely  equivalent  with  starch  (Gijzen  1985).  Conversion  of  starch  to  glucose 
increases  the  dry  matter  with  a  factor  10/9,  whereas  the  transport  of  glucose 
costs  dry  matter,  i.e.,  36/38.  Thus,  the  total  transport  “cost”  of  downward  trans¬ 
location  is  a  factor  CVT  =  1.05  (10/9  x  36/38).  Measured  data  on  translocation 
are  extremely  scarce  for  terrestrial  plants  and  absent  for  aquatic  plants.  Trans¬ 
location  proved  to  be  29  percent  of  net  production  in  cassava  (Gijzen  1985)  and 
35  percent  in  certain  potato  varieties  (Kooman  1996).  Estimates  of  translocation 
in  submersed  plants  vary  from  19  percent  of  net  production  in  sea  grasses 
(Wetzel  and  Heckles  1996)  to  about  40  percent  in  hydrilla  (Best  and  Boyd  1996). 

In  MILFO,  TRANS  follows  a  hyperbolic  relationship  initially  set  to  35  per¬ 
cent  (TRAFAC)  of  net  production  by  the  senescing  plant  cohort,  multiplied  by 
CVT,  and  decreasing  exponentially  to  zero  with  concomitantly  decreasing  bio¬ 
mass  of  the  translocating  plant  cohort  and  increasing  biomass  of  the  successive 
growing  plant  cohort. 

Translocation  is  described  by  the  following  equation: 

TRANSl  =  CVT  X  GPHOT  x  ((TWLG2+TWSG2+TWRG2)I 

(TWLG1+TWSG1+TWRG1+TWLG2+TWSG2+TWRG2))  x 
TRAFAC 
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Figure  3.  Relational  diagram  illustrating  photosynthesis,  respiration,  and 
biomass  formation  in  milfoil 
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where 


TRANSl  =  translocation  rate  cohort  1  (g  CH2O  m'^  d"*) 

CVT  =  conversion/transport  factor  (relative) 

GPHOT  =  daily  total  gross  CHjO  assimilation  rate  of  community 
(gCH^Om-^d-i) 

TWLGl  or  2  =  total  weight  of  green  leaves  cohort  1  or  cohort  2  (g  DW  m'^ 

TWSGl  or  2  =  total  weight  of  green  stems  cohort  1  or  cohort  2  (g  DW  m‘^ 

TWRGl  0x2  =  total  weight  of  live  roots  cohort  1  or  cohort  2  (g  DW  m'^ 
TRAFAC  =  translocation  factor  (relative) 

Senescence  refers  to  the  loss  of  capacity  to  carry  out  essential  physiological 
processes  and  to  the  loss  of  biomass.  The  fundamental  processes  involve 
physiological  aging  and  protein  (enzyme)  breakdown.  These  processes  are 
difficult  to  quantify.  It  is  known  that  hormones  are  important  messengers  in  this 
context,  but  not  how  they  precisely  act.  High  temperature  usually  accelerates 
senescence. 

In  MILFO,  a  mechanistic  approach  to  senescence  has  been  chosen  by  setting 
the  death  rate  at  a  certain  fraction  of  plant  biomass  lost  per  day  once  the  condi¬ 
tions  for  growth  deteriorate.  The  timing  and  value  of  relative  death  rate  (RDR) 
of  plant  cohorts  1  and  2  have  been  derived  from  field  observations  on  shoot  bio¬ 
mass  in  Lake  Wingra,  Wisconsin  (Adams  and  McCracken  1974).  For  plant 
cohort  3,  timing  and  relative  death  rate  of  the  plant  cohorts  1  and  2  performed 
well,  and  length  of  the  third  senescence  period  turned  out  to  be  similar  to  that  of 
plant  cohort  2. 

The  timing  of  onset  of  senescence  was  found  by  running  the  model  repeatedly 
with  different  development  rates,  base  and  reference  temperatures.  Thus,  initia¬ 
tion  of  senescence  for  cohort  1  was  set  at  DVS  1.631,  for  cohort  2  at  DVS  2.001, 
and  for  cohort  3  at  3.501. 

The  value  for  the  relative  death  rate  of  the  plants  was  found  by  applying  the 
same  differential  equation  as  commonly  used  for  simple  exponential  growth,  to 
describe  exponential  decrease  in  biomass  after  flowering,  with  a  negative  speci¬ 
fic  decrease  rate  (Thomley  and  Johnson  1990b;  Hunt  1982).  An  RDR  of 
0.042  d’’  was  calculated  following  this  approach. 

The  value  for  the  relative  death  rate  of  the  rhizome/root  crown  system  was  set 
at  0.00042  d'*  as  described  in  the  section  Wintering  and  Sprouting  of  Rhizomes/ 
Root  Crowns  and  Growth  of  Sprouts  to  Water  Surface. 
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A  relational  diagram  illustrating  translocation  and  senescence  is  shown  in 
Figure  4. 


Figure  4.  Relational  diagram  illustrating  translocation  and  senescence  following 
anthesis  in  milfoil 
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Choice  of  Parameter  Values 


A  relatively  simple  simulation  model  like  MILFO  includes  parameter  values  that  can 
be  defined  with  varying  certainty.  Most  parameters  have  been  calculated/estimated  from 
published  literature  (Table  3).  Only  development  rate  in  relation  to  3  °C  day-degree 
sum  and  base  temperature  have  been  calibrated  by  running  the  model.  The  choice  of 
parameter  values  has  been  detailed  in  the  preceding  sections  of  this  chapter. 
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Table  3 

Parameter  Values  Used  in  MILFO 


Parameter 


Abbreviation  Value 


Morphology,  Phenological  Cycle,  and  Development 


Reference 


First  Julian  day  number 


Development  rate  as  function  of  temperature 


Base  temperature  for  juvenile  plant  growth 


DAYEM 


DVR(T)* 


TEASE 


0.015-0.022 


C 


Calibrated 


Calibrated 


Maximum  Biomass  and  Plant  Density 


Maximum  biomass 

AMIN1  TGW 

2,283  gDWnr^ 

4 

Plant  density 

NPL 

11  m'^ 

4 

Wintering  and  Sprouting  of  RhIzomes/Root  Crowns  and  Growth  of  Sprouts  to  Water  Surface 

Critical  rhizome  weight 

CRRIZ 

12gDWm-^ 

8 

Initial  rhizome  weight 

IWGRIZ 

50gDWm-“ 

10 

Relation  coefficient  rhizome/root  crown  weight-stem  length 

RCSHST 

12  m  g  DW'^ 

3.17 

Relative  death  rate  of  rhizomes 

RDRIZ 

0.00042  d’ 

18 

Relative  conversion  rate  of  rhizome/root  crown  weight  into  plant  material 

ROC 

0.0576  g  CHjO. 
g  DW’ d  ' 

3 

I  Light  and  Photosynthesis  of  Piants 

Daytime  temperature  effect  on  AMX  as  function  of  DVS 

AMTMPfT) 

0-1 

13,  11 

Potential  COg  assimilation  rate  at  light  saturation  for  shoot  tips 

AMX 

0.0165  gCOj. 
g  DW’  It’ 

2,16 

Conversion  factor  for  translocated  dry  matter  into  CHgO 

CVT 

1.05 

9 

Initial  light-use  efficiency  for  shoot  tips 

EE 

0.000011  gCOj  J  ’ 

9 

Reflection  coefficient  of  Irradiance  at  water  surface 

RC 

0.06 

5 

Reduction  factor  to  relate  AMX  to  water  pH 

REDAM 

0.5 

7,15 

Reduction  factor  for  AMX  to  account  for  senescence  plant  parts  over 
vertical  vegetation  axis 

REDF(T) 

1.0 

User  def. 

Plant  species  specific  light-extinction  coefficient 

K{T) 

0.006  g  DW’ 

13 

Water  type  specific  light-extinction  coefficient 

LfD 

1.15-2.00  m  ’ 

7 

Thickness  per  plant  layer 

TL 

0.1  m 

14 

Maintenance,  Growth,  and  Assimilate  Partitioning  of  Plants 

Dry  matter  allocation  to  each  plant  layer 

DMPC(T) 

0-1 

1 

Leaf  dry  matter  allocation  to  each  plant  layer 

FLCT) 

0.50 

1 

_ 

Fraction  of  total  dry  matter  increase  allocated  to  leaves 

FLV(D 

0.47 

1 

Fraction  of  total  dry  matter  increase  allocated  to  roots 

FRT(T) 

0.06 

1 

Fraction  of  total  dry  matter  increase  allocated  to  stems 

FST{D 

0.47 

1 

Factor  accounting  for  effect  of  dally  effective  temperature  on  maintenance 
respiration 

TEFF(T) 

0-12 

12 

(Continued) 


28 


Chapter  3  Model  Processes 


liable  3  (Concluded)  ll 

Parameter 

Abbreviation 

Value 

Reference 

Flowering,  Translocation,  and  Senescence 

Relative  death  rate  of  leaves  (on  DW  basis) 

RDR(T) 

0.042  d-' 

2 

Relative  death  rate  of  stems  and  roots  (on  DW  basis) 

RDS(T) 

0.042  d' 

2 

Translocation  (part  of  net  photosynthetic  rate) 

TRAFAC 

0.35 

5 

Site  Information 

Lag  period  between  water  and  air  temperature 

DELAY 

7d 

User  def. 

Water  depth  (=  rooting  depth) 

DEPTH 

1.5  m 

User  def. 

Total  live  dry  weight  measured  (field  site) 

TGWM{T) 

-,  g  DM  m-" 

User  def. 

Daily  water  temperature  (field  site) 

WTMPfT) 

User  def. 

Harvesting 

Harvesting 

HAR 

Oor  1 

User  def. 

Harvesting  day  number 

HARDAY 

1-365 

User  def. 

Harvesting  depth  (measured  from  water  surface  in  m) 

HARDER 

0.1  m<DEPTH 

User  def. 

Notes:  1.  Adams,  Titus,  and  McCracken  1974;  2.  Adams  and  McCracken  1974;  3.  Bowes,  Holaday,  and  Haller  1979;  4.  Budd, 
Lillie,  and  Rasmussen  1995;  5.  Golterman  1975;  6.  Kooman  1995;  7.  Lee  and  Kluesener  1972;  8.  Madsen  1997;  9.  Penning  de 
Vries  and  Van  Laar  1 982a,  b;  1 0,  Smith  and  Adamsi  986;  1 1 .  Stanley  and  Nailor  1 972;  12.  Thornley  and  Johnson  1 990a;  1 3. 

Titus  and  Adams  1979a;  14.  Titus  et  al.  1975;  15.  Titus  and  Stone  1982;  16.  Van,  Haller,  and  Bowes  1976;  17.  Van  der  Zweerde 

1 981 ;  1 8.  Voqt,  Vogt,  and  Bloomfield  1 991 .  *,  Calibration  function. 
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4  Performance  Tests 


Simulated  and  Measured  Behavior  of  a  Milfoil 
Community  in  Lake  Wingra,  Wisconsin 

Nominal  run 

The  seasonal  changes  in  biomass  of  plant  shoots  and  roots  and  of  the 
rhizome/root  crown  system  as  simulated  by  MILFO  are  shown  in  Figure  5A  and 
B.  Simulated  shoot  biomass  compared  well  with  shoot  biomass  measured  in 
Lake  Wingra  (Figure  5C).  Peak  biomass  appeared  to  be  reached  somewhat 
earlier  in  the  simulation  than  found  in  the  lake;  however,  the  latter  may  be  due  to 
the  low  frequency  of  field  observations  (no  measurements  between  September 
and  November).  The  simulated  biomass  of  the  rhizome/root  crown  system 
showed  two  maxima  per  year.  Variation  was  within  the  range  found  in  a  milfoil 
community  in  the  same  lake  in  later  years  (Smith  and  Adams  1986). 

Simulated  transport  of  carbohydrates  was  substantial  in  the  beginning  and  at 
the  end  of  the  growing  period  of  each  plant  cohort,  when  carbohydrate  remobili¬ 
zation  from  the  rhizome/root  crown  system  supports  growth  of  the  sprouts,  and 
carbohydrate  translocation  from  plant  organs  supports  the  filling  of  the  rhizome/ 
root  crown  system,  respectively  (Figure  6).  Carbohydrate  transport  could  be  in 
the  same  range  as  net  assimilation  at  the  beginning  and  end  of  the  growth  season 
(Figure  7).  Maintenance  respiration  was  usually  considerably  lower  than  assimi¬ 
lation  as  well  as  carbohydrate  transport,  as  can  be  expected  at  the  relatively  low 
water  temperatures  (Figure  7). 

Running  the  model  with  the  low  assimilate  requirements  suggested  to  be  typi¬ 
cal  for  submersed  aquatic  vegetation  by  Spencer,  Ryan,  and  Ksander  (1997) 
showed  that  peak  biomass  of  milfoil  shoots  would  increase  by  a  factor  of  2, 
oscillations  in  biomass  of  the  rhizome/root  crown  system  would  be  greater,  and 
final  biomass  of  the  rhizome/root  crown  system  would  be  increased  (Figure  8). 
However,  as  indicated  in  Chapter  3,  the  opinion  of  these  authors  is  that  a  con¬ 
struction  cost  of  0.99  to  1.11  for  milfoil  plant  tissues  is  on  the  low  side,  taking 
the  usually  high  N  concentrations  of  shoots  into  consideration. 
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43“  08'  N).  Measured  shoot  biomass  data  1970  (Adams  and  McCracken 
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Figure  6.  Simulated  behavior  of  carbohydrate  flow  through  plant  compartments 
(Carbohydrate  remobilization  and  upward  transport  from  rhizomes/ 
root  crowns  are  used  for  initial  growth  of  each  plant  cohort.  Down- 
Ward  carbohydrate  translocation  into  rhizomes/root  crowns  occurs 
during  anthesis  and  senescence  of  each  plant  cohort  (Initial  biomass 
and  climatological  data  as  in  nominal  run)) 


Running  the  model  for  the  same  year  and  lake,  but  with  only  rhizomes/root 
crowns  initially  present  (Figure  9B),  showed  that  plant  biomass  of  both  cohorts 
was  greatly  reduced  and  critical  weight  of  the  rhizome/root  crown  system  was 
reached  more  often  than  with  initial  plant  biomass  present  (Figure  9A;  nominal 
run).  This  large  difference  in  peak  biomass  is  due  to  the  inability  of  the  plant 
community  to  fully  capture  the  high  spring  irradiance  at  this  latitude  of  43  °N 
without  wintering  shoots.  Thus,  wintering  shoots  provide  a  distinct  advantage 
for  this  plant  species. 

Running  the  model  with  (24-hr  average)  air  temperatures  lagging  7  days 
behind  water  temperatures  or  measured  water  temperatures  as  forcing  variables 
yielded  similar  biomass  values,  despite  the  fact  that  instantaneous  assimilation 
rates  varied  less  with  water  temperatures  than  with  air  temperatures,  and  assimi¬ 
lation  rates  had  shifted  somewhat  in  time  (Figure  10).  This  illustrates  the 


Chapter  4  Performance  Tests 


d> 

0) 

^  c. 

CO  o 

§  2  -b 

«  S'  <v‘ 
S  E 

^  i 

CD 


CD 

%  g 

'5 


d 

3? 

O 

.o? 


10 


8  H 


6H 


4H 


2  H 


S;  it; 


Simulated  values^ 
Lake  Wingra,  Wl 


-  net  assimilation 


5:. 


}  :-:• 


*  <»  • 


maintenance 


0 


yf^'v 


maintenance 

rhizomes 


100 


200 


Time  (d) 


T 

300 


400 


Figure  7.  Simulated  rates  of  daily  net  assimilation  and  maintenance  respiration 
of  a  milfoil  community  in  Lake  Wingra,  Wisconsin  (Initial  biomass  and 
climatological  data  as  in  nominal  run) 


usefulness  of  inclusion  of  both  temperature  options  in  the  model,  facilitating  use 
of  the  model  by  users  who  do  not  possess  a  full  data  set  of  water  temperatures 
for  the  water  body  for  which  they  desire  to  run  the  model.  It  has  to  be  cautioned, 
however,  that  the  relationship  between  the  temperatures  of  air  and  of  each  water 
body  concerned  may  differ  since  temperatures  within  each  water  body  are  influ¬ 
enced  by  catchment  morphometry,  wind  speed,  fetch,  mixing  processes,  and 
upward  seepage,  etc.  In  the  experience  of  these  authors,  however,  a  lag  period  of 
7  days  between  air  and  water  temperatures  usually  described  this  relationship 
well  for  shallow  water  bodies  (up  to  5-6  m),  without  large  inflows  of 
groundwater. 


Effects  of  year-to-year  differences  in  ciimate 

The  model  was  run  with  initial  biomass  values  and  local  climatological  data 
as  inputs  for  a  different  year,  1972  (Figure  11).  A  run  with  water  temperatures 
of  a  previous  (1970)  year  yielded  less  biomass  (Figure  llA)  than  actually 
measured  (Figure  IIC).  A  run  with  air  temperatures  of  1972,  in  contrast,  yielded 
less  biomass  for  the  first  plant  cohort,  but  similar  biomass  as  measured  for  the 
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as  in 
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Figure  9.  Simulated  biomass  with  (A)  and  without  (B)  initial  plant  biomass  and  measured  (C)  shoot  biomass  of  a  milfoil  community  in  Lake 
Wingra,  Wisconsin  ((A)  Initial  biomass  values  and  climatological  data  as  in  nominal  run;  (B)  initial  plant  biomass  zero  and 
rhizomes/root  crowns  as  in  nominal  run.  Measured  shoot  biomass  1970  (Adams  and  McCracken  1974).  Climatological  data  as  in 
nominal  run) 


Figure  10.  Simulated  photosynthetic  rates  of  a  milfoil  community  in  Lake 

Wingra,  Wisconsin,  with  water  or  air  temperatures  as  input  (Initial 
biomass  and  climatological  data  as  in  nominal  run) 


second  cohort  (Figure  IIB).  However,  irradiances  in  both  years,  1970  and  1972, 
differed  in  that  total  irradiance  and,  consequently,  temperature  sum  were  higher 
in  1970  than  in  1972  (particularly  in  spring),  and,  thus,  higher  biomass  produc¬ 
tion  was  to  be  expected  using  water  temperatures  of  1970.  This  leads  one  to 
believe  that  the  early  peak  biomass  value  measured  in  1972  is  an  overestimate. 
The  latter  suggestion  is  supported  by  the  fact  that  the  measured  biomass  level 
could  neither  be  attained  by  running  the  model  with  considerably  decreased 
light-extinction  coefficients  tentatively  indicative  for  the  clear  water  phase  in 
spring,  which  is  typical  for  this  lake. 

Simulated  and  Measured  Behavior 

of  a  Milfoil  Community  at  Other  Latitudes 

To  investigate  whether  the  model  was  able  to  simulate  behavior  of  a  milfoil 
community  at  other  sites,  runs  were  made  for  a  site  at  a  more  southern  latitude. 
Lake  Guntersville,  Alabama.  Behavior  of  milfoil  in  this  lake  is  particularly 
interesting  because  the  lake  is  long,  oriented  from  north-east  to  south-west,  and 
situated  at  a  latitude  around  34  °N,  being  very  close  to  tropical  (33  °N).  Bio¬ 
mass  of  milfoil  communities  in  this  lake  has  been  described  as  having  a  high 
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variation  in  time  and  space  (Grace  and  Wetzel  1978;  Stanley  et  al.  1976);  unfor¬ 
tunately,  in  these  descriptions,  no  attention  was  paid  to  differences  in  latitude  of 
the  various  sites  within  the  lake  nor  to  local  differences  in  temperature  or  other 
environmental  factors.  It  was  mentioned,  however,  that  flowering  and  the  sub¬ 
sequent  sloughing  period  were  less  predictable  in  southern  locales  than  in  nor¬ 
thern  ones. 

Initial  plant  biomass  values  measured  at  a  site  in  this  lake  studied  in  1990 
were  very  low,  possibly  because  of  grass  carp  herbivory  the  previous  year,^  and 
initial  rhizome/root  crown  mass  has  been  set  at  the  critical  value  of  12  g  DW  m"^. 
Rooting  depth  in  the  simulation  was  kept  at  1.5  m,  although  in  reality  water 
depth  may  have  varied  over  1.0  +  0.7  m  within  a  year  (Stanley  et  al.  1976). 

Simulated  biomass  of  the  first  plant  cohort  remained  low.  Only  one  apparent 
biomass  peak  could  be  distinguished,  which  originated  from  the  second  plant 
cohort.  Simulated  shoot  biomass  coincided  in  timing  with  measured  shoot  bio¬ 
mass,  but  the  simulated  peak  was  lower  than  the  measured  one  (Figure  12).  The 
latter  difference  may  be  a  consequence  of  temporal  decreases  in  water  depth  dur¬ 
ing  the  year;  depth  was  kept  constant  in  the  simulation,  leading  to  an  underesti¬ 
mate  of  simulated  plant  biomass.  Relatively  small  changes  in  water  depth  can 
cause  large  changes  in  net  assimilation  and  biomass  production  (See  Chapter  6). 

To  investigate  which  consequences  a  warm  year  for  the  milfoil  community  in 
this  lake  might  have,  when  three  instead  of  two  plant  cohorts  are  expected  to  be 
active,  a  model  run  was  made  with  the  same  initial  biomass  and  climatological 
data  and  a  third  cohort  active  (Figure  13).  It  turned  out  that  in  one  year  far 
higher  shoot  biomass  values  of  approximately  950  g  DW  m'^  could  be  generatedj 
similar  in  timing  and  value  to  maximum  biomass  values  reported  for  the  nearby 
Melton  Hill  Lake  (Stanley  et  al.  1976),  with  rhizome/root  crown  biomass  accum¬ 
ulating  towards  the  end  of  the  year.  However,  similar  biomass  values  could  also 
be  reached  earlier  in  the  year  when  higher  initial  (nominal)  biomass  values  were 
used  as  input,  and  only  two  cohorts  active;  in  the  latter  case,  biomass  peaks  of 
both  plant  cohorts  appeared,  and  rhizome/root  crown  biomass  was  well  above 
the  critical  level  but  not  accumulating. 

Simulation  of  a  milfoil  community  in  the  Kashmir  lakes,  India,  demonstrated 
that  only  one  maximum  in  shoot  biomass  was  generated  (Figure  14),  with  a  value 
somewhat  higher  than  the  measured  range  of  288  to  640  g  DW  m‘^  and  a 
rhizome/root  crown  biomass  varying  over  a  range  close  to  the  measured  range  of 
32  to  160  (Zutschi  and  Vass  1973).  The  higher  simulated  shoot  biomass  may  be 
due  to  the  use  of  climatological  data  from  Patancheru,  which  is  located  more 
south,  and  thus  warmer,  than  the  Kashmir  lakes  from  which  the  measured 
biomass  values  originated  (Patancheru  17  °N,  Kashmir  lake  32  °N);  however, 
more  northern  climatological  data  from  India  were  not  available. 


‘  Personal  Communication,  M.S.  Stewart,  U.S.  Army  Engineer  Waterways  Experiment  Station, 
Vicksburg,  MS. 
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Figure  13.  Simulated  biomass  of  shoots,  roots,  and  rhizomes/root  crowns  of  a  milfoil  community  in  Lake  Guntersville,  Alabama,  composed  of 
two  or  three  plant  cohorts  and  using  various  initial  biomass  values  (Measured  initial  biomass  vlaues — cf.  M.  S.  Stewart 
unpublished  (1990)  and  nominal  ones  as  for  nominal  run.  Climatological  data  1990,  Lake  Guntersville,  Alabama) 


Simulated  plant  biomass  Simulated  rhizomes  1  Aboveoround  plant  biomass 

^‘atancheru,  Jndia  Patancheru,  India  Patancneru,  India 
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lakes,  India  (longitude  74  to  80°  E,  latitude  32°  N;  Initial  values  on  biomass  of  shoots  and  rhizomes/root  crowns  mean  of  range 
published  for  1970s  (Zutschi  and  Vass  1973).  Climatological  data,  1978,  Patancheru,  India  (longitude  78°  28'  E,  latitude  17°  27'  N). 
Measured  shoot  biomass  mean  of  range  published  for  1970s  (Zutschi  and  Vass  1973)) 


It  was  investigated  whether  milfoil  benefits  from  adaptation  to  the  tropics  by 
producing  thinner  leaves.  This  was  done  because  a  higher  leaf-surface  area:dry 
weight  ratio  (K-value)  has  been  found  for  milfoil  in  Japan  (0.01  m^  g  DW'‘; 
Ikusima  1970)  than  in  Wisconsin  (0.006  m^  g  DW‘^;  Titus  and  Adams  1979a).  It 
turned  out  that  timing  was  very  similar  and  simulated  plant  biomass  about 
10  percent  higher  using  the  higher  K-value  (data  not  shown). 

Running  the  model  with  nominal  biomass  values  and  climatological  data  typi¬ 
cal  for  sites  representative  for  temperate,  temperate  to  tropical,  and  tropical  cli¬ 
mates  (Figure  15)  indicated  that  (a)  in  all  climates  one  clear  biomass  peak  is 
generated;  (b)  only  in  a  temperate  climate  the  biomass  peak  of  both  first  and 
second  cohorts  can  be  distinguished;  that  is,  from  biomass  values  alone;  flower¬ 
ing  coinciding  with  every  biomass  maximum  is  always  a  suitable  indicator,  but  it 
is  often  not  noted  in  biomass  studies;  (c)  peak  biomass  is  expected  to  be  highest 
in  the  tropics;  that  highest  biomass  values  have  been  found  at  northern  latitudes 
may  be  because  most  biomass  studies  on  aquatic  plants  have  been  performed  at 
the  latter  latitudes  and  biomass  data  from  tropical  areas  are  extremely  scarce; 
and  (d)  end-of-year  accumulation  of  rhizome/root  crown  biomass  usually  occurs 
in  tropical,  but  not  in  temperate  climates;  that  is,  when  three  plant  cohorts  are 
active. 
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5  Sensitivity  Anaiysis 


A  sensitivity  analysis  of  a  simulation  model  is  required  to  assess  the 
parameters  likely  to  strongly  affect  model  behavior.  The  current  analysis  was 
based  on  the  effect  of  a  change  in  a  parameter  when  all  other  parameters  are  kept 
the  same.  As  reference  level,  the  nominal  parameter  values  were  chosen  as 
presented  in  Table  3,  under  Lake  Wingra,  Wisconsin,  conditions  at  1.5-m  water 
depth.  In  a  1-year  simulation  starting  with  50  g  DW  m‘^  biomass  for  both  plants 
and  rhizome/root  crown  system,  the  value  of  the  parameter  under  study  was 
changed  (Table  4).  The  results  were  compared  with  those  of  a  nominal  run. 

Each  parameter  was  once  increased  by  20  percent  and  once  decreased  by 
20  percent.  The  relative  sensitivity  (RS)  of  a  parameter  was  then  defined  as  the 
relative  change  in  the  variable  on  which  the  effect  was  tested  divided  by  the 
relative  change  in  the  parameter  (Ng  and  Loomis  1984).  The  effects  of 
10  parameters  on  two  variables,  representing  plant  biomass  aspects,  were  tested. 
A  model  variable  is  considered  sensitive  to  a  change  in  the  value  of  a  parameter 
at  RS  >  0.5  and  <-0.5.  The  current  sensitivity  analysis  was  performed  over  a 
1-year  period. 

(yield-  -  yield I  yield ^ 

RS  —  - 

(param-  -  param^l param^ 


where 

yieldi  =  value  at  parameter  value  i 
yields  =  value  at  reference  parameter  value 
paranii  and  param^  as  above 

Maximum  plant  biomass  proved  most  sensitive  to  changes  in  potential  CO2 
assimilation  at  light  saturation  for  shoot  tips  and  very  sensitive  to  changes  in 
light-use  efficiency.  This  is  not  surprising  because  the  model  is  based  on  carbon 
flow  through  the  plant.  Changes  in  plant  density  did  affect  maximum  plant 
biomass,  but  far  less  than  photosynthetic  activity  at  light  saturation  and  light-use 
efficiency.  Most  parameter  changes,  except  in  critical  rhizome  weight, 
influencing  rhizome/root  crown  biomass  affected  maximum  plant  biomass 


44 


Chapter  5  Sensitivity  Anaiysis 


substantially,  for  example,  initial  rhizome  weight,  conversion  rate  into  plant 
material,  and  translocation  rate.  In  general,  the  same  parameters  as  those  for 
maximum  plant  biomass  were  important  determinants  of  end-of-year 
rhizome/root  crown  biomass,  with  potential  CO2  assimilation  at  light  saturation, 
light-use  efficiency,  and  relative  death  rate  exhibiting  the  largest  effects.  This 
illustrates  the  utmost  importance  of  the  rhizome/root  crown  system  for  local 
survival  and  biomass  production  of  milfoil. 

Earlier  or  later  flowering  biotypes  are  suited  to  different  environments.  The 
effect  of  flowering  date  can  be  tested  with  the  model  by  varying  the  development 
rate  of  the  vegetation.  Slower  rates  represent  later  and  faster,  earlier  biotypes. 
Development  rate  slower  or  faster  than  the  nominal  rate  leads  to  lower  biomass. 
Faster  development  leads  to  a  shorter  growing  season  and  less  vegetative  dry 
matter,  incomplete  light  interception,  and  lower  carbohydrate  availability  for 
organ  formation.  At  the  same  time,  however,  the  rate  of  organ  formation 
increases,  but  the  duration  of  each  organ  formation  shortens.  Intuitive  prediction 
of  biotype  behavior  under  such  highly  variable  climatic  conditions  is  therefore 
hazardous.  The  model  shows  some  promise  in  being  able  to  reproduce  some  of 
these  complex  responses  of  the  vegetation  and  may  be  useful  in  evaluating  long¬ 
term  implications  of  differences  in  development  rate. 

Although  as  far  as  is  known,  no  publications  exist  on  what  the  temperature 
requirements  of  aquatic  plants  are  to  traverse  development  from  anthesis  to 
senesced  state;  differences  in  postanthesis  development  rates  for  several  wheat 
and  rice  cultivars  are  known  to  be  small  and  have  little  effect  on  yield  (V an 
Keulen  1976). 

Maximum  plant  biomass  proved  only  sensitive  to  a  decrease  in  preanthesis 
development  rate,  while  end-of-year  rhizome/root  crown  biomass  was  sensitive 
to  any  change  in  preanthesis  or  postanthesis  development  rate. 
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Table  4 

Relative  Sensitivity  of  Two  Modei  Variabies  to  Deviations  in  Parameter  Vaiues  from 

Their  Nominai  Vaiues  (As  presented  in  Tabie  3)  (Resuits  were  obtained  in  a  1-year 
simuiation  under  Lake  Wingra,  Wisconsin,  1970  conditions,  starting  with  both  piant  and 
rhizome  biomass  being  50  g  DW  m'^ 

Relative  Sensitivity 

Parameter  Name 

Parameter  Value 

Maximum  Live  Plant 

Biomass 

End-of-Year 
Rhizome/Root  Crown 
Biomass 

Potential  CO2  assimilation  rate  at  light  saturation 
for  shoot  tips 

0.0165 

0.0200 

1.96 

2.00 

0.0149 

1.97 

2.02 

Light-use  efficiency 

0.000011 

0.000013 

1.10 

1.14 

0.000008 

1.22 

1.25 

Relative  death  rate  leaves,  stems,  and  roots 

0.042 

0.050 

-0.62 

-1.01 

0.034 

-0.77 

-1.36 

Initial  rhizome  weight 

50 

60 

0.20 

0.17 

40 

0.22 

0.18 

Critical  rhizome  weight 

12 

14.4 

0 

0.05 

9.6 

0 

0.06 

Relative  conversion  rate  of  rhizomes  Into  plant 
material 

0.0576 

0.069 

0.19 

0.17 

0.046 

0.21 

0.18 

Translocation  rate 

0.35 

0.42 

-0.13 

0.57 

0.28 

-0.14 

|0.72 

Plant  density 

'11 

13 

0.16 

-0.79 

9 

-0.16 

0.79 

Preanthesis  development  rate 

0.015 

0.018 

i  -0.23 

-0.81 

0.012 

-0.26 

-0.99 

Postanthesis  development  rate 

0.015 

0.018 

-0.69 

-0.89 

0.012 

-0.79 

0.66 
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6  Environmental  Factor 
Analysis 


The  impacts  of  various  changes  in  environmental  factors  were  assessed  using 
the  relative  sensitivity  of  the  affected  variables  as  “measure.”  For  this  purpose, 
parameter  changes  were  based  on  value  ranges  taken  from  literature,  which 
sometimes  differed  more  than  20  percent  from  the  nominal  parameter  value 
given  in  Table  3. 


Climate 

Climate  greatly  affects  plant-species  distribution,  phenological  cycle,  and 
biomass  production.  MILFO  can  be  used  to  calculate  climate  change  effects  on 
the  chronological  timing  of  the  phenological  events  and  on  biomass  production. 
It  cannot  be  used  to  assess  climate  change  effects  on  (a)  plant-species 
distribution  and  (b)  the  phenological  cycle  itself  since  the  phenological  cycle  has 
been  used  for  calibration  (see  Chapter  3).  Running  the  model  under  more 
southern  climatological  conditions,  i.e.,  changing  the  latitude  from  43  to  34°  N 
demonstrated  that  end-of-year  rhizome/root  crown  biomass  is  far  more  sensitive 
to  this  climate  change  than  maximum  plant  biomass  (Table  5). 


Light-Reflection  Coefficient  at  Water  Surface 

The  irradiance  reflected  at  the  water  surface  usually  averages  about  6  percent 
daily.  The  values  of  this  parameter  tested  were  0  and  1.  Reflection  may 
theoretically  have  the  value  0  when  no  reflection  occurs  at  a  90  °  incoming  angle 
of  the  radiation  on  a  completely  calm  water  surface  (wind  and  wave  action  are 
minimal).  The  highest  value  of  1  may  occur  at  a  close  to  180  °  incoming  angle 
of  the  radiation  and  at  very  rough  water  surfaces. 

Increasing  the  light  reflection  coefficient  to  1  annihilated  plant  biomass 
within  the  year.  That  nevertheless  low  RS  values  were  found  (Table  5)  is  an 
artifact  of  the  calculation  method  employed.  Decreasing  the  light-reflection 
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Table  5 

Environmental  Factor  Analysis,  Expressed  as  Relative  Sensitivity  of  Two  Model 
Variables  to  Deviations  in  Parameter  Values  from  Their  Nominal  Values  (As  presented 
in  Table  3)  (Results  were  obtained  in  a  1-year  simulation  under  Lake  Wingra, 
Wisconsin,  1970  conditions,  starting  with  both  plant  and  rhizome/root  crown  biomass 
being  50  g  DW  m'^ ) _ 


Parameter  Name 

Parameter  Value 

Relative  Sensitivity 

Maximum  Live  Plant 
Biomass 

End-of-Year  Rhizome/Root 
Crown  Biomass 

Climate 

Uke  Wingra  (1980) 

Latitude  43°  N 

- 

- 

Lake  Guntersville  (1990) 

Latitude  34°  N 

-0.01 

0.89 

Light-reflection  coefficient  at  water  surface 

0.06 

1.00  (+1567%) 

-0.05 

-0.04 

0.00*  (-100%) 

-0.07 

-0.07 

Light-extinction  coefficient  water  column 

1.80 

2.16  (+20%) 

-0.92 

-1.01 

1.44  (-20%) 

-1.12 

-1.01 

Water  depth 

1.5 

1.8  (+20%) 

-0.31 

-0.31 

1.2  (-20%) 

-0.34 

-0.33 

[Note:  To  enable  calculation  of  the  RS,  a  very  low  value  of  0.000001  was  used. 


coefficient  greatly  increased  maximum  biomass  and  end-of-year  rhizome/root 
crown  biomass  (Table  5). 

Light-Extinction  Coefficient  of  Water  Coiumn 

A  light-extinction  coefficient  of  on  average  1.80  m'^  is  used  for  nominal  runs 
of  the  model  (Lake  Wingra,  Wisconsin). 

Changing  the  light-extinction  coefficient  of  the  water  column  demonstrated 
large  effects  on  maximum  plant  and  end-of-year  rhizome/root  crown  biomass. 
The  often  relatively  small  effect  of  an  increase  in  light-extinction  coefficient 
relative  to  the  nominal  value  may  be  due  to  (a)  the  high  nominal  value  and 
(b)  the  spatial  distribution  of  milfoil  plant  biomass  with  typically  61  percent  in 
the  upper  0.5-m  water  layer.  A  nominal  value  of  2  m'^  has  been  found  typical  for 
eutrophic  fen  lakes  where  submersed  vegetation  can  just  persist  (Best,  De  Vries, 
and  Reins  1985).  The  large  effect  of  a  decrease  in  light-extinction  coefficient 
can  largely  be  explained  by  greatly  increased  growth  of  the  first  plant  cohort. 
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boosting  the  rhizome/root  crown  system  by  translocation,  thus  providing  a  better 
start  for  the  subsequent  plant  cohort(s)  and  resulting  in  a  higher  peak  biomass. 


Water  Depth 

MILFO  has  been  calibrated  for  a  water  depth  of  1.5  m,  the  rooting  depth  of 
an  extensively  studied  milfoil  community  in  Lake  Wingra,  Wisconsin.  The 
model  has  the  capability  to  respond  to  fluctuations  in  water  level  with  year,  by 
(re)distributing  plant  biomass  over  the  desired  water  depth  (number  of  water 
layers;  see  Chapter  3).  This  technique  for  biomass  distribution  over  the  vertical 
axis  of  the  community  works  well  and  gives  realistic  outcomes  over  a  depth 
range  of  0.5  to  6  m. 

Running  MILFO  at  an  increased  or  decreased  water  depth  showed  similar, 
relatively  small  effects  on  maximum  plant  and  end-of-year  rhizome/root  crown 
biomass  (Table  5). 

The  RS  of  peak  plant  biomass  and  end-of-year  rhizome/root  crown  biomass 
to  changes  in  water  depth  was  relatively  small  and  far  lower  than  to  changes  in 
the  light-extinction  coefficient. 

The  current  sensitivity  and  environmental  analyses  give  indications  of  the 
sensitivity  of  maximum  plant  biomass  and  end-of-year  rhizome/root  crown  sys¬ 
tems  for  variations  in  plant  parameters  and  environment  over  a  1-year  period.  It 
is  to  be  expected,  however,  that  the  small  changes  that  occurred  over  this  rela¬ 
tively  short  period  will  increase  with  time  and  that  these  extrapolations  in  time 
will  yield  information  on  the  likelihood  for  plant  populations  to  ultimately  per¬ 
sist  or  become  extinct.  Particularly,  increases  in  water  turbidity  because  of 
increased  phytoplankton  or  periphyton  growth  stimulated  by  eutrophication, 
increased  erosion/resuspension,  and  seasonal  herbivory  have  been  mentioned  as 
decisive  for  the  persistence  of  submersed  plant  populations. 
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7  Application  Possibilities 


MILFO  can  be  used  to  assess  behavior  of  a  milfoil  community  under  various 
site-specific  and  climatological  conditions  as  demonstrated  in  Chapters  4, 5,  and 
6,  and  it  can  be  run  with  user-specified  input  values  for  plant  and  rhizome/root 
crown  biomass. 

Effects  of  man-made  control  activities  like  harvesting  at  different  times  and  at 
various  water  depths  can  be  calculated  also  (Table  6).  Thus,  in  the  latter  case  it 
can  be  used  as  a  tool  for  aquatic  plant  management  agencies.  From  this  table  it 
can  be  concluded  that  harvesting  at  the  end  of  May  to  a  water  depth  of  0.8  m 
requires  removal  of  a  relatively  low  amount  of  biomass,  but  yields  the  lowest 
peak  biomass  and  end-of-year  rhizome/root  crown  biomass.  This  situation  can 
be  seen  as  favorable  to  control  milfoil.  In  contrast,  harvesting  later  in  the  year 
requires  removal  of  relatively  more  plant  biomass  or  allows  for  a  relatively 
higher  end-of-year  rhizome/root  crown  biomass.  Removing  only  the  top  layer  of 
the  plant  community  later  in  the  year  may  even  lead  to  increased  maximum  plant 
and  end-of-year  rhizome/root  crown  biomass,  probably  because  of  a  temporarily 
higher  light  penetration  within  the  community. 


Table  6 

Effects  of  Mechanical  Harvesting  Date  and  Depth  on  Plant  and  Rhizome/Root  Crown 
Biomass  (Results  were  obtained  in  a  1-year  simulation  under  Lake  Wingra,  Wisconsin, 
1970  conditions,  starting  with  both  plant  and  rhizome/root  crown  biomass  being  50  g 

DW  m*) 

Harvest  Time 

Harvest 

Depth 

m 

Live  Piant 
Biomass 

28  August 
g  DW  m-" 

Preharvest 
Biomass 
g  DW  m  " 

Postharvest 
Biomass 
g  DW  m'* 

Day  with 
Zero  Piant 
Biomass 

End-of-Year  Rhizome/ 
Root  Crown  Biomass 
g  DW  m^ 

End  of  May 

0.8 

84 

100 

21 

>365 

13 

End  of  June 

0.8 

143 

106 

5 

>365 

18 

End  of  July 

0.8 

49 

171 

8 

>365 

14 

End  of  July 

0.1 

300 

171 

141 

>365 

41 

End  of  August 

0.8 

258 

244 

11 

>365 

18 

End  of  September 

0.8 

258 

100 

4 

>365 

30 

End  of  October 

0.8 

258 

33 

1 

>365 

33 

50 


Chapter  7  Application  Possibilities 


The  current  version  of  MILFO  has  been  developed  as  a  stand-alone  simula¬ 
tion  model.  It  can  be  relatively  easily  modified  to  communicate  with  ecosystem 
models  because  it  is  written  in  FORTRAN??  and  its  structure  is  simple.  It  is 
planned  to  link  MILFO  to  a  Geographical  Information  System  through  an  appro¬ 
priate  interface  like  AEGIS+  (Luyten  et  al.  1994).  To  facilitate  use  of  the  cur¬ 
rent  model,  a  user  manual  has  been  prepared  (Best  and  Boyd,  in  preparation). 


Chapter  7  Application  Possibilities 
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8  Discussion 


The  current  model  gives  a  reasonable  description  of  the  dynamics  in  plant 
and  rhizome/root  crown  biomass  of  an  established  milfoil  population  under  a 
variety  of  field  conditions.  As  can  be  expected,  the  model  is  very  sensitive  to 
environmental  changes  affecting  the  light  climate  and,  consequently,  the  carbon 
flow  through  the  plant. 

Extinction  of  light  by  periphyton  has  not  been  included  in  MILFO  because 
(a)  the  plant  canopy  tends  to  be  at  the  water  surface  during  most  of  the  growth 
season,  (b)  irradiance  in  the  euphotic  zone  of  the  plant  canopy  (upper  layers)  is 
often  saturating  (i.e.,  >600  uE  cm‘^  s'^;  Van,  Haller,  and  Bowes  1976),  and  (c)  no 
field  data  on  periphyton  biomass  concomitant  with  photosynthetic  activity  are 
available  at  this  time.  Light  attenuation  by  periphyton  is  expected  to  have  large 
effects  on  submersed  macrophyte  species  with  most  of  their  biomass  concen¬ 
trated  just  above  the  hydro-soil  (like  Ceratophyllum  demersum;  Best  and  Dassen 
1987;  Best  and  Jacobs  1990)  and  macrophytes  with  biomass  that  usually  remains 
below  the  water  surface  (like  Vallisneria  americana;  Titus  and  Adams  1979a). 

Senescence,  resulting  in  decreasing  photosynthetic  activity  in  aging  plant 
parts,  has  been  included  into  the  model  formulation.  Although  data  quantifying 
these  effects  in  milfoil  were  available  (Adams  and  McCracken  1974;  Adams, 
Titus,  and  McCracken  1974),  running  the  model  demonstrated  that  virtually  no 
effect  on  peak  biomass  was  noticeable,  probably  largely  because  of  the  typical 
umbrella-type  biomass  distribution  over  the  water  column,  with  not  only  most 
biomass  in  the  upper  portion  of  the  community  but  also  most  young  plant  parts. 
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♦ _ * 

*  SUBROUTINE  MODEL  * 

*  Authors:  Elly  Best  &  Will  Boyd  * 

*  Date  :  16  July  1997  * 

*  * 


*  FORMAL  PARAMETERS:  (I=input,0=output,C=control,IN=init,T=time) 


*  name 

type 

meaning 

units 

class  * 
_ ♦ 

*DELT 

R4 

Time  step  of  integration 

d 

I  * 

*DOY 

R4 

Day  number  within  year  of  simulation  (REAL) 

d 

I  * 

♦FILEIN 

C* 

Name  of  file  with  input  model  data 

- 

I  ♦ 

♦FINTIM 

R4 

Finish  time  of  simulation  (=day  number) 

d 

I  * 

*IDOY 

14 

Day  number  witliin  year  of  simulation  (INTEGER) 

d 

I  ♦ 

*  ITASK 

14 

Task  that  subroutine  should  perform 

- 

I  * 

♦lUNITD 

14 

Unit  of  input  file  with  model  data 

- 

I  ♦ 

♦lUNITO 

14 

Unit  of  output  file 

- 

I  * 

♦lUNITL 

14 

Unit  number  for  log  file  messages 

- 

I  * 

♦lYEAR 

14 

Year  of  simulation  (INTEGER) 

y 

I  * 

*LAT 

R4 

Latitude  of  site 

dec.degr. 

I  * 

♦LONG 

R4 

Longitude  of  site 

dec.degr. 

I  ♦ 

♦ELEV 

R4 

Elevation  of  site 

m 

I  ♦ 

♦  OUTPUT 

L4 

Flag  to  indicate  if  output  should  be  done 

I  * 

♦RAIN 

R4 

Daily  amount  of  rainfall 

mm/d 

I  ♦ 

♦RDD 

R4 

Daily  shortwave  radiation 

J/m2/d 

I  * 

♦STTIME 

R4 

Start  time  of  simulation  (=day  number) 

d 

I  ♦ 

♦TERMNL 

L4 

Flag  to  indicate  if  simulation  is  to  stop 

- 

I/O* 

♦  TMMN 

R4 

Daily  minimum  temperature 

degrees  C 

I  * 

♦TMMX 

R4 

Daily  maximum  temperature 

degrees  C 

I  * 

♦VP 

R4 

Early  morning  vapour  pressure 

kPa 

I  ♦ 

♦  WN 

R4 

Daily  average  wind  speed 

m/s 

I  * 

♦  WSTAT 

C6 

Status  code  from  weather  system 

- 

I  * 

♦WTRTER 

L4 

Flag  whether  weather  can  be  used  by  model 

- 

0  * 

♦YEAR 

R4 

Year  of  simulation  (REAL) 

y 

I  * 

♦  * 

*  Fatal  error  checks:  if  one  of  the  characters  of  WSTAT  =  '4',  indicates  missing  weather  * 

*  Warnings  :  none  * 

*  Subprograms  called:  models  as  specified  by  the  user  * 

*  File  usage  :  lUNTTO.IUNITD+LIUNITO.IUNITO+LIUNITL  * 

*  _ * 


SUBROUTINE  MODEL  (ITASK ,  lUNITD,  lUNTTO,  lUNITL,  FBLEIN, 
&  OUTPUT,  TERMNL, 

&  DOY  ,  IDOY  ,  YEAR  ,  lYEAR, 

&  TIME  ,  STTIME,  FINTIM,  DELT , 

&  LAT  ,  LONG  ,  ELEV  ,  WSTAT ,  WTRTER, 

&  RDD  ,  TMMN  ,  TMMX  ,  VP  ,  WN,  RAIN) 

*  Title  of  the  program 

*  <Fill  in  your  title  here> 

*  MILFO 

IMPLICIT  REAL  (A-Z) 

* — -Formal  parameters 
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INTEGER  ITASK ,  lUNITD,  lUNITO,  lUNTTL,  IDOY,  lYEAR 

LOGICAL  OUTPUT,  TERMNL,  WTRTER 

CHARACTER*(*)  FILEIN,  WSTAT 

REAL  DOY,  YEAR,  TIME,  STTIME,  FINTIM,  DELT 

REAL  LAT,  RDD,  TMMN,  TMMX,  VP,  WN,  RAIN 

REAL  TMAX(365),  TMIN(365) 


* - Standard  local  declarations 

INTEGER  IWVAR,ITOLD,IDAY,  DDELAY 
CHARACTER  WUSED*6 


—State  variables,  initial  values  and  rates 
REALDVS  ,NUL  ,DVR  ,TMPSUM,TMP2 
REAL  TWLGl  ,TWSG1  .TWRGl  ,IWLG1  ,IWSG1,IWRG1 
REALTWLDl  ,TWSD1  ,TWRD1  ,IWLD1  ,IWSD1,IWRD1 
REAL  TWLG2  ,TWSG2  ,TWRG2  ,IWLG2  ,IWSG2,IWRG2 
REAL  TWLD2  ,TWSD2  ,TWRD2  ,IWLD2  ,IWSD2,IWRD2 
REALTWLG3  ,TWSG3  ,TWRG3  ,IWLG3  ,IWSG3,IWRG3 
REAL  TWLD3  ,TWSD3  ,TWRD3  ,IWLD3  ,IWSD3,IWRD3 
REAL  NGLVl  ,NGST1  ,NGRT1  ,DLV1  ,DST1  ,DRT1 
REAL  NGLV2  ,NGST2  ,NGRT2  ,DLV2  ,DST2  ,DRT2 
REAL  NGLV3  ,NGST3  ,NGRT3  ,DLV3  ,DST3  ,DRT3 
REAL  TWLVG  ,TWLVD  ,TWSTG  ,TWSTD  ,TWRTG,TWRTD 
REAL  TWGRIZ,TWRIZD,IWGRIZ,IWRIZD,TGRIZ 


* — ^Model  parameters 

REALAMX  ,CVT  ,DAYEM , DELAY  , DEPTH , EE 
REALHAR  ,HARDAY,HARDEP 
REAL  NPL  ,RC  ,RCSHST,REDAM ,  RDRIZ 
REAL  ROC  ,TBASE  ,TRAFAC,TL 


—Auxiliary  variables 

REAL  AMAX  ,  AlvITMP ,  ASRQ  ,  COSLD  ,  WTMP 
REAL DAVTMP, DAY  ,DAYL  ,DDTMP,DSO 
REAL  DSINB  ,  DSINBE,  DTEFF  ,  DTGA  ,  FGROS 
REALFLV  ,FRT  ,FRT1  ,FRT2  ,FST 
REALGLV  ,GPHOT,GRT  ,  GST  ,  GTW 
REALMAINT,MAINTS,PI  ,  RDR  ,RDS 
REAL  REMOBl,  REMOB2,  REMOB3,  SC  ,  SINLD 
REAL  SUM  ,  TEFF  ,  TGW  ,  TGWM  ,  TRANSl 
REAL  TRANS2,  TRANS3,  TREMOB,  TW  ,  WLV 
REALWST  ,WRT  ,  MAINRT,  YRNUM 


* - AFGEN  functions 

REALAMTMPT 
INTEGER  IMAMTM,  ILAMTM 
PARAMETER  (1MAMTM=  40) 
DIMENSION  AMTMPT(IMAMTM) 
REALFLT 

INTEGER  IMFLT,  ILFLT 
PARAMETER  (IMFLT  =  40) 
DIMENSION  FLT  (IMFLT) 
REALFLVT 
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INTEGER  IMFLVT,  ILFLVT 
PARAMETER  (IMFLVT=  40) 
DIMENSION  FLVT  (IMFLVT) 
REALFRTT 

INTEGER  IMFRTT,  ILFRTT 
PARAMETER  (IMFRTT=  40) 
DIMENSION  FRTT  (IMFRTT) 
REALFSTT 

INTEGER  IMFSTT,  ILFSTT 
PARAMETER  (IMFSTT=  40) 
DIMENSION  FSTT  (IMFSTT) 
REAL  LT,  KT 
INTEGER  IMN1,ILT,IKT 
PARAMETER  (IMNl  =  730) 
DIMENSION  LT(IMNI),  KT(IMNl) 
REALTGWMT 
INTEGER  IMMEAS,  ILMEAS 
PARAMETER  (IMMEAS  =  40) 
DIMENSION  TGWMT(IMMEAS) 
REAL  RDRT 

INTEGER  IMRDRT,  ILRDRT 
PARAMETER  (IMRDRT=  40) 
DIMENSION  RDRT  (IMRDRT) 
REALRDST 

INTEGER  IMRDST,  ILRDST 
PARAMETER  (IMRDST  =  40) 
DIMENSION  RDST  (IMRDST) 
REALTEFTT 

INTEGER  IMTEFF,  ILTEFF 
PARAMETER  (IMTEFF=  40) 
DIMENSION  TEFFT  (IMTEFF) 
REALWTMPT 

INTEGER  IMWTMP,  ILWTMP 
PARAMETER  (IMWTMP  =  730) 
DIMENSION  WTMPT  (IMWTMP) 


*  - ^Used  functions 

REAL  LINT  ,INSW 
SAVE 

DATA  ITOLD  /4/ 

*  - Code  for  the  use  of  HDD,  TMMN,  TMMX,  VP,  WN,  RAIN  (in  that  order) 

*  A  letter  'U'  indicates  that  the  variable  is  used  in  calculations 
DATA  WUSED/’UUU— 7 


* - Check  weather  data  availability 

IF  aTASK.EQ,1.0R.n’ASK.EQ.2.0R.ITASK.EQ.4)  THEN 
DO  10  IWVAR=1,6 

* — “Is  there  an  error  in  the  IWVAR-th  weather  variable  ? 

IF  (WUSED(IWVAR:IWVAR).EQ.’LP  .AND. 

&  WSTAT(IWVAR:IWVAR).EQ.'4')  THEN 
WTRTER  =  .TRUE. 

TERMNL  =  .TRUE. 

ITOLD  =  ITASK 
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RETURN 
END  IF 
10  CONTINUE 
END  IF 

IF  (ITASK.EQ.l)  THEN 

* 

*  INITIALIZATION  SECTION 

♦  ======================== 

* - Send  title  to  output  file 


* - Open  input  file 

CALL  RDINIT  (lUNITD,  lUNITL,  FILEIN) 


* - ^Read  1st  value  in  MODEL.DAT  file ...  year  number 

CALL  RDSREA  ('YRNUM  '.YRNUM ) 


* - ^Read  initial  states 

CALL  RDSREA  (‘CRRIZ  ',CRRIZ ) 
CALL  RDSREA  (TWGRIZ'JWGRIZ) 
CALL  RDSREA  (TWRIZD'JWRIZD) 
CALL  RDSREA  (TWLD 1  '.IWLD 1 ) 
CALL  RDSREA  (TWLD2  MWLD2  ) 
CALL  RDSREA  (TWLD3  MWLD3  ) 
CALL  RDSREA  (TWLGl  MWLGl ) 
CALL  RDSREA  (TWLG2  ■,IWLG2  ) 
CALL  RDSREA  (TWLG3  ',IWLG3  ) 
CALL  RDSREA  (TWRDl  '.IWRDl ) 
CALL  RDSREA  (TWRD2  ',IWRD2 ) 
CALL  RDSREA  (TWRD3  MWRD3  ) 
CALL  RDSREA  (TWRGl  MWRGl ) 
CALL  RDSREA  (TWRG2  MWRG2  ) 
CALL  RDSREA  (TWRG3  ’,IWRG3  ) 
CALL  RDSREA  (TWSDl  ',IWSD1 ) 
CALL  RDSREA  ('IWSD2  ’JWSD2  ) 
CALL  RDSREA  (TWSD3  MWSD3  ) 
CALL  RDSREA  ('IWSGl  MWSGl ) 
CALL  RDSREA  (TWSG2  ',IWSG2  ) 
CALL  RDSREA  (TWSG3  ',IWSG3  ) 
CALL  RDSREA  ('NUL  ',NUL  ) 
CALL  RDSREA  ('REMOBl',REMOBl) 


♦ - Read  model  parameters 

CALL  RDSREA  ('AMX  ',AMX  ) 
CALL  RDSREA  ('CVT  '.CVT  ) 

CALL  RDSREA  ('DAYEM  '.DAYEM ) 
CALL  RDSREA  ('DELAY  '.DELAY  ) 
CALL  RDSREA  ('DEPTH  '.DEPTH  ) 
CALL  RDSREA  ('EE  '.EE  ) 

CALL  RDSREA  CHAR  '.HAR  ) 

CALL  RDSREA  (HARDAY'.HARDAY) 
CALL  RDSREA  (’HARDEP'.HARDEP) 
CALL  RDSREA  ('NPL  '.NPL  ) 

CALL  RDSREA  ('RC  '.RC  ) 
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CALL  RDSREA  ('RCSHST’,RCSHST) 
CALL  RDSREA  ('REDAM  ',REDAM  ) 
CALL  RDSREA  ('RDRIZ  ',RDRIZ ) 
CALL  RDSREA  ('ROC  ',ROC  ) 
CALL  RDSREA  ('XBASE  '.XBASE ) 
CALL  RDSREA  ('XL  ',XL  ) 

CALL  RDSREA  ('XRAFAC'.XRAFAC) 


* - ^Read  AFGEN  functions 

CALL  RDAREA  ('AMXMPX',AMXMPX,IMAMXM,ILAMXM) 
CALL  RD AREA  CFLX  ',FLX  ,IMFLX  ,ILFLX  ) 

CALL  RDAREA  ('FLVX  '.FLVX  ,IMFLVX,n.FLVX) 

CALL  RDAREA  ('FSXX  '.FSXX  ,IMFSXX,ILFSXX) 

CALL  RDAREA  ('FRXX  '.FRXX  ,IMFRXX,ILFRXX) 

CALL  RDAREA  CKX  ',KX  ,IMN1  ,IKX  ) 

CALL  RDAREA  ('LX  ',LX  ,IMN1  ,ILX  ) 

CALL  RDAREA  ('RDRX  ',RDRX  ,IMRDRX,ILRDRX) 

CALL  RDAREA  ('RDSX  ',RDSX  ,IMRDSX,ILRDSX) 

CALL  RDAREA  CXEFFX  '.XEFFX  ,IMXEFF,E.XEFF) 

CALL  RDAREA  ('XGWMX  '.XGWMX  ,IMMEAS,ILMEAS) 
CALL  RDAREA  ('WXMPX  '.WXMPX  ,IM\VXMP,ILWXMP) 

***  INIXIAL  CALCULAXIONS  *** 


*  - Initially  known  variables  to  output 

*  Send  title(s)  to  OUXCOM 

*  - Initialize  state  variables 

*  Start  at  tlie  beginning  of  tlie  developmental  cycle 

DVS  =NUL 
XMPSUM  =  Nl]L 


+ - ^Delay  variable  set  from  a  REAL  to  an  INXEGER 

DDELAY  =  DELAY 


* - ^Initialize  weights  of  plant  organs 

IF  (YRNUM  .EQ.  l.)XHEN 


XWLDl 

XWLD2 

XWLD3 

XWLGl 

XWLG2 

XWLG3 

XWSDl 

XWSD2 

XWSD3 

XWSGl 

XWSG2 

XWSG3 

XWRDl 

XWRD2 

XWRD3 

XWRGl 


=  IWLD1 
=  IWLD2 
=  IWLD3 
=  IWLG1 
=  IWLG2 
=  I\VLG3 
=  IWSDl 
=  IWSD2 
=  IWSD3 
=  IWSGl 
=  IWSG2 
=  IWSG3 
=  IWRD1 
=  1WRD2 
=  IWRD3 
=  IWRG1 
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TWRG2  =IWRG2 
TWRG3  =IWRG3 
TWGRIZ  =  IWGRIZ 
TWRIZD  =  IWRIZD 
ENDIF 

ELSE  IF  (ITASK.EQ.2)  THEN 
***  RATES  OF  CHANGE 


* - Weights  of  plant  organs 

WLV  =  TWLGl  +  TWLDl  +  TWLG2  +  TWLD2  +  TWLG3  +  TWLD3 

WST  =  TWSGl  +  TWSDl  +  TWSG2  +  TWSD2  +  TWSG3  +  TWSD3 

WRT  =  TWRGl  +  TWRDl  +  TWRG2  +  TWRD2  +  TWRG3  +  TWRD3 

TGW  =  TWLG1+TWSG1+TWRG1+TWLG2+TWSG2+TWRG2+TWLG3+TWSG3+TWRG3 


*  - ^Total  live  weight  never  >2283  g  DW  /  m2;  cf.  Budd  et  al. 

TGW  =  AMINl  (TGW,  2283.) 

*  - ^Initialize  rhizome  weight 

TWRIZ  =  TWGRIZ  -  TWRIZD 

***  RATE  CALCULATIONS  ♦** 


*  - ^Julian  day  number 

DAY  =  l.+MOD  (TIME-1.,365.) 

*  - If  water  temperatures  are  available,  temperature  dependent  processes  are  related  to  water 

*  temperature;  otherwise  they  are  related  to  air  temperature  with  a  lag  period  in  day(s)  to  be  chosen 

*  by  substituting  number  given  for  DELAY  in  MODEL.DAT 

WTMP  =  LINT  (WTMPT,ILWTMP,DAY) 

IDAY  =  DAY 

TMAX(IDAY)  =  TMMX 
TMIN(IDAY)  =  TMMN 
IF  (DAY  .LE.  DDELAY)  THEN 

DAVTMP  =  0.5  *  (TMAX(1)+TMIN(1)) 

DDTMP  =  TMAX(l)  -  0.25  *  (TMAX(l)-TMIN(l)) 

ELSE 

DAVTMP  =  0.5  *  (TMAXaDAY-DDELAY)+TMIN(IDAY-DDELAY)) 

DDTMP  =  TMAX(IDAY-DDELAY)  -  0.25  * 

&  (TMAX(IDAY-DDELAY)-TMIN(IDAY-DDELAY)) 

ENDIF 

IF  (DAVTMP  .LT.  5.0)DAVTMP  =  5.0 

IF  (WTMP  .GT.  0.0)  THEN 
DAVTMP- WTMP 
DDTMP  =  WTMP 
ENDIF 
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TEFF  =  LINT(TEFFT,ILTEFF,DDTMP) 


— ^Measured  total  live  plant  dry  weight 

TGWM  =  LINT  (TGWMT,ILMEAS,DAY) 

- Call  to  SBRT  ASTRO  to  introduce  day  length  into  MAIN  for  tentative 

relationship  REMOBl-DAYL;  otherwise  tliis  call  can  be  made  at  *ASTRO 
CALL  ASTRO 

$(DAY,LAT,SC,DSO,SINLD,COSLD,DAYL,DSINB,DSINBE) 

— Calculation  of  dry  matter  and  its  partitioning  over  tlie  plant  organs 
TW  =  TGW+(TWLD1+TWSD1+TWRD1+TWLD2+TWSD2+TWRD2+TWLD3+TWSD3+TWRD3) 
FLV  =  LINT(FLVT  ,ILFLVT,DVS) 

FST  =  LINT(FSTT  ,ILFSTT,DVS) 

FRT  =  LINT(FRTT  .ILFRTT,DVS) 

FL  =  LINT(FLT  ,ILFLT  ,DVS) 


“Growth  of  plant  organs,  maintenance  respiration  and  translocation 
Calculation  assimilate  requirement  for  plant  organ  formation. 
ASRQ  =1.54 

MAINTS  =  0.016*TWLG1+0.01*TWSG1+0.015*TWRG1 
MAINT  =MAINTS*TEFF 
MAINRT  =  0.005  *  TWGRIZ  *  TEFF 


— Carbohydrate  behavior;  remobilization  from  rhizomes  for  plant  formation  at  proper  day  length  and 
temperature  conditions  (presently  related  to  DVS);  translocation  from  above-ground  biomass 
(i.e. ’plants')  to  rhizomes,  provided  plants  are  present 
TRANS  1  =  0.0 
TREMOB  =  0.0 

TWGRIZ  =  AMAXl  (CRRIZ,  TWGRIZ) 

TGRIZ  =  TWGRIZ 

TWRIZD  =  INTGRL  (TWRIZD,  RDRIZ ,  BELT) 

IF  (DVS.GE.0.376  .AND.  DVS.LT.1.0)  THEN 
IF  (GPHOT  .LT.  MAINT)THEN 
IF  (TWGRIZ  .GT.  CRRIZ)  THEN 
TREMOB  =  INTGRL  (TREMOB,  REMOBl,  DELT) 

REMOBl  =  ROC*  TWGRIZ 
ELSE 

WRITE  (*,*)  'Vegetation  is  dying' 

REMOBl  =  0.000001 
ENDIF 

ENDIF 

ELSE 

REMOBl  =  0.000001 
ENDIF 


—Relative  death  rates 

RDR  =  INSW  (DVS-1.001,0.,LINT  (RDRT,ILRDRT,DDTMP)) 
RDS  =  INSW  (DVS-1.001,0.,LINT  (RDST,ILRDST,DDTMP)) 
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*  - ^Development  rates 

IF(DAVTMP  XT.  3.0)  THEN 
DVR  =  0.0 

ELSE  IF  (DVS.LE.l.)  THEN 
DVR  =  0.022*DAVTMP/30 
ELSE  IF  (DVS.LE.6.0)  THEN 
DVR  =  O.OI5*DAVT1VIP/30 
ENDIF 

*  - Calculation  of  effective  daytime  temperature 

DTEFF  =  AMAX1(0.,DAVTMP-TBASE) 


Calculation  of  dead  plant  material  per  organ 
DLVl  =  TWLGl  *  RDR 
DSTl  =  TWSGl  *  RDR 
DRTl  =  TWRGl  *  RDR 

■Shoot  photosynthesis  at  tight  saturation,  and  daytime  temperature  effect  on  shoot  photosynthesis 
AMAX  =  AMAX1(0.00001,AMX  *  AMTMP) 

AMAX  =  AMAX*  REDAM 

AMTMP  =  LINT(AMTMPT,ILAMTM,DDTMP) 


* - ^Before  calling  TOTASS,  determine  light  extinction  coefficients  of  plant  material  (K)  and  water  (L) 

L  =  LINT(LT,ILT,TIME) 

K  =  LINT(KT,IKT,DVS) 


* - Daily  total  gross  assimilation 

CALL  TOTASS 

$  (SC,DAYL,SINLD,COSLD,DSINBE,RDD,RC,L,K,AMAX,EE, 
$  TL,DEPTH,RCSHST,TGW,FGROS,FL,WLV,WST, 

$  DAY,HAR,HARDAY,HARDEP,DTGA,IRS) 


*  - If  harvesting  takes  place,  weights  various  plant  organs  must  be  recalculated 

*  (TWLG1,TWSG1,TWRG1,TW) 

IF(HAR  .EQ.  1.  AND.  DAY  .EQ.  HARD  AY)  THEN 
TWLGl  =FLV*TGW 
TWSGl  =  FST  *  TGW 
TWRGl  =  FRT  *  TGW 

TW  =  TGW+(TWLDl+TWSDl+TWRDl+TWLD2+TWSD2+TWRD2+TWLD3-t-TWSD3+TWRD3 
ENDIF 

*  - Conversion  assimilated  C02  to  CH20 

GPHOT  =  DTGA  *  30./44. 

*  - Induction  of  flowering  at  DVS=1;  flowering  occurs  10  days  after  induction 

*  Induction  of  flowering,  translocation  and  senescence  occur  simultaneously 
IF  (DVS  .GE.  1.0)THEN 


*  - If  there  is  no  above-ground  plant  biomass  present,  TRANS  1  must  stay  at  zero;  otherwise  it  gets 

*  a  value 

IF  ((TWLGl+TWSGl+TWRGl)  .GT.  0.0)THEN 
IF  (TGW  .GE.  TWGRIZ)  THEN 

TRANSl  =  CVT  *  GPH0T*((TWLG1+TWSG1+TWRG1)/ 
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$(TWLG1+TWSG1+TWRG1+TWLG2+TWSG2+TWRG2))  *  TRAFAC 

ENDIF 

ENDIF 

*  - If  tliere  is  no  plant  biomass  present,  REMOB  1  must  stay  at  zero 

IF  ((TWLGl+TWSGl+TWRGl)  .EQ.  0.0)REMOB1  =  0.0 
ENDIF 

*  The  total  weight  of  the  green  rhizomes  increases  by  translocation  (TRANSl  from  CHORTl)  and 

*  decreases  by  remobilization  and  rhizome  maintenance  respiration 
TWGRIZ  =  INTGRL  (TWGRIZ,  -(((-TRANSl  +REMOB1+  MAINRT)/1.242) 

$  +(RDRIZ*TGRIZ)),  BELT) 

*  - Total  and  net  growtli  rates 

GTW  =  ((REMOB  1*CVT)  +  GPHOT  -  TRANSl  -  MAINT)/ASRQ 
GRT  =  FRT*GTW 
GST  =  FST  *  GTW 
GLV  =  FLV*GTW 

NGLVl  =  GLV  -  DLVl 
NGSTl  =  GST  -  DSTl 
NGRTl  =  GRT-DRT1 

CALLCHORT2  (AMAX  .AMTMP  ,AMTMPT,AMX  ,ASRQ  ,COSLD  ,CRRIZ , 

&  CVT  ,DAVTMP,DAY  ,DAYL  ,DDTMP,DELT  ,DEPTH, 

&  DLV2  ,DRT2  ,DSINB  ,DSINBE,DST2  ,DTEFF  ,DTGA  , 

&  DVR  ,DVS  ,EE  ,FGROS,FL  ,FLT  ,FLV  , 

&  FLVT  ,FRT  ,FRTT  ,FST  ,FSTT  ,GLV  ,GPHOT, 

&  GRT  ,GST  ,GTW  ,HAR  ,HARDAY,HARDEP,IKT  , 

&  ILAMTM,ILFLT  ,ILFLVT,ILFRTT,ILFSTT,ILRDRT,ILRDST, 

&  ILT  ,K  ,KT  ,L  ,LAT  ,LT  , 

&  MAINRT,MAINT  ,NGLV2  ,NGRT2  ,NGST2  ,RC  , 

&  RCSHST,RDD  ,RDR  ,RDRIZ,RDRT  ,RDST  , RED  AM, 

&  REMOB l,REMOB2,ROC  ,SC  ,SINLD  ,TBASE  ,TEFF  , 

&  TGRIZ,TGW  ,TIME  ,TL  ,TRAFAC,TRANS1,TRANS2, 

&  TREMOB,TW  ,TWLD1  ,TWLD2  ,TWLG1  ,TWLG2  ,TWRD1 , 

&  TWRD2  ,TWRG1,  TWRG2  ,TWGRIZ,TWRIZD,TWSD1 , 

&  TWSD2  ,TWSG1,  TWSG2  ,WLV  ,WST  ,IRS) 

IF  (LAT  .LE.  33)  CALL  CHORT3 

&  (AMAX , AMTMP  ,AMTMPT,AMX  ,ASRQ  ,COSLD  ,CRRIZ , 

&  CVT  ,DAVTMP,DAY  ,DAYL  ,DDTMP,DELT  ,DEPTH, 

&  DLV3  ,DRT3  ,DSINB  ,DSINBE,DST3  ,DTEFF  ,DTGA  , 

&  DVR  ,DVS  ,EE  ,FGROS,FL  ,FLT  ,FLV  , 

&  FLVT  ,FRT  ,FRTT  ,FST  ,FSTT  ,GLV  .GPHOT, 

&  GRT  ,GST  ,GTW  ,HAR  ,HARDAY,HARDEP,IKT  , 

&  ILAMTM,ILFLT  ,ILFLVT,ILFRTT,ILFSTT,ILRDRT,ILRDST, 

&  ILT  ,K  ,KT  ,L  ,LAT  ,LT  , 

&  MAINRT,MAINT  ,NGLV3  ,NGRT3  ,NGST3  ,RC  , 

&  RCSHST,RDD  ,RDR  ,RDRIZ  ,RDRT  ,RDST  , REDAM , 

&  REMOBl,REMOB3,ROC  ,SC  ,SINLD  ,TBASE  ,TEFF  , 

&  TGRIZ.TGW  .TIME  ,TL  ,TRAFAC,TRANS1,TRANS2, 

&  TRANS3,TREMOB,TW  .TWLDl  ,TWLD2  ,TWLD3  .TWLGl , 
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TWLG2  ,TWLG3  ,TWRDl  ,TWRD2  ,TWRD3  ,TWRG1  ,TWRG2 , 
TWRG3  ,TWGRIZ,TWRIZD,TWSD1  ,TWSD2  ,TWSD3  , 

TWSGl  .TWSG2  ,TWSG3  ,WLV  ,WST  ,IRS) 


* - ^Finish  conditions 

IF(DVS  .GT.  6.0  .OR.  DAY  .EQ.  365.)  TERMNL  =  .TRUE. 


* - Output  section 

IF  (OUTPUT)  THEN 

CALL  OUTDAT  (2,0,'DAVTMP',DAVTMP) 
CALL  OUTDAT  (2,0, DAYL  ',DAYL  ) 
CALL  OUTDAT  (2,0,'DDTMP  ',DDTMP ) 
CALL  OUTDAT  (2,0,DTEFF  •,DTEFF ) 
CALL  OUTDAT  (2,0,’DTGA  ',DTGA  ) 
CALL  OUTDAT  (2,0,DVS  •,DVS  ) 

CALL  OUTDAT  (2,0,'FGROS  ',FGROS  ) 
CALL  OUTDAT  (2,0,'GPHOT  '.GPHOT ) 
CALL  OUTDAT  (2,0,'IRS  ',IRS  ) 

CALL  OUTDAT  (2,0,'MAINRr,MAINRT) 
CALL  OUTDAT  (2,0,'MAINT  '.MAINT ) 
CALL  OUTDAT  (2,0,'REMOB1',REMOB1) 
CALL  OUTDAT  (2,0,’REMOB2',REMOB2) 
CALL  OUTDAT  (2,0,'REMOB3‘,REMOB3) 
CALL  OUTDAT  (2,0,’TEFF  ',TEFF  ) 

CALL  OUTDAT  (2,0,'TGW  ',TGW  ) 
CALL  OUTDAT  (2,0,'TGWM  ',TGWM  ) 
CALL  OUTDAT  (2,0,’TMPSUM',TMPSUM) 
CALL  OUTDAT  (2,0,TRANS1',TRANS1) 
CALL  OUTDAT  (2,0,'TRANS2',TRANS2) 
CALL  OUTDAT  (2,0,'TRANS3',TRANS3) 
CALL  OUTDAT  (2,0,'TREMOB',TREMOB) 
CALL  OUTDAT  (2,0,TW  ',TW  ) 

CALL  OUTDAT  (2,0,TWGRIZ',TWGRIZ) 
CALL  OUTDAT  (2,0,'TWLDl  ’,TWLD1 ) 
CALL  OUTDAT  (2,0,TWLD2  ',TWLD2 ) 
CALL  OUTDAT  (2,0,TWLD3  ',TWLD3  ) 
CALL  OUTDAT  (2,0,TWLG1  ',TWLG1 ) 
CALL  OUTDAT  (2,0,'TWLG2  ',TWLG2  ) 
CALL  OUTDAT  (2,0,'TWLG3  ■,TWLG3  ) 
CALL  OUTDAT  (2,0,TWLVD  ',TWLVD  ) 
CALL  OUTDAT  (2,0,TWLVG  ’,TWLVG ) 
CALL  OUTDAT  (2,0,’TWRDl  ■,TWRD1 ) 
CALL  OUTDAT  (2,0,TWRD2  ■,TWRD2  ) 
CALL  OUTDAT  (2,0,'TWRD3  •,TWRD3  ) 
CALL  OUTDAT  (2,0,'TWRGl  ’,TWRG1  ) 
CALL  OUTDAT  (2,0,'TWRG2  •,TWRG2  ) 
CALL  OUTDAT  (2,0,TWRG3  ',TWRG3  ) 
CALL  OUTDAT  (2,0,'TWRIZ  ’,TWRIZ ) 
CALL  OUTDAT  (2,0,TWRIZD',TWRIZD) 
CALL  OUTDAT  (2,0,'TWRTG  •,TWRTG  ) 
CALL  OUTDAT  (2,0,'TWSDl  ',TWSD1 ) 
CALL  OUTDAT  (2,0,'TWSD2  ',TWSD2  ) 
CALL  OUTDAT  (2,0,TWSD3  ',TWSD3  ) 
CALL  OUTDAT  (2,0,'TWSGl  ■,TWSG1 ) 


Appendix  A  Model  Listing 


A11 


CALL  OUTDAT  (2,0,TWSG2  ',TWSG2  ) 
CALL  OUTDAT  (2,0,TWSG3  ’,TWSG3  ) 
CALL  OUTDAT  (2,0,'TWSTD  ■,TWSTD ) 
CALL  OUTDAT  (2,0,TWSTG  ’,TWSTG  ) 
CALL  OUTDAT  (2,0,’WTMP  ’,WTMP  ) 

END  IF 

ELSE  IF  (ITASK.EQ.3)  THEN 
INTEGRATION 


DVS  =INTGRL(DVS  ,DVR  ,DELT) 
TMPSUM  =  INTGRL  (TMPSUM,DTEEF  ,DELT) 
TWLDI  =  INTGRL  (TWLDI  ,DLV1  ,DELT) 
TWLD2  =  INTGRL  (TWLD2  ,DLV2  ,DELT) 
TWLD3  =  INTGRL  (TWLD3  ,DLV3  ,DELT) 
TWLGI  =  INTGRL  (TWLGl  ,NGLV1  ,DELT) 
TWLG2  =  INTGRL  (TWLG2  ,NGLV2  ,DELT) 
TWLG3  =  INTGRL  (TWLG3  ,NGLV3  ,DELT) 
TWLGI  =  AMAXl  (0.0,  TWLGl) 

TWLG2  =  AMAXl  (0.0,  TWLG2) 

TWLG3  =  AMAXl  (0.0,  TWLG3) 

TWSDl  =  INTGRL  (TWSDl  ,DST1  ,DELT) 
TWSD2  =  INTGRL  (TWSD2  ,DST2  ,DELT) 
TWSD3  =  INTGRL  (TWSD3  ,DST3  ,DELT) 
TWSGl  =  INTGRL  (TWSGl  ,NGST1  ,DELT) 
TWSG2  =  INTGRL  (TWSG2  ,NGST2  ,DELT) 
TWSG3  =  INTGRL  (TWSG3  ,NGST3  ,DELT) 
TWSGl  =  AMAXl  (0.0,  TWSGl) 

TWSG2  =  AMAXl  (0.0,  TWSG2) 

TWSG3  =  AMAXl  (0.0,  TWSG3) 

TWRGl  =  INTGRL  (TWRGl  ,NGRT1  ,DELT) 
TWRG2  =  INTGRL  (TWRG2  ,NGRT2  ,DELT) 
TWRG3  =  INTGRL  (TWRG3  ,NGRT3  ,DELT) 
TWRGl  =  AMAXl  (0.0,  TWRGl) 

TWRG2  =  AMAXl  (0.0,  TWRG2) 

TWRG3  =  AMAXl  (0.0,  TWRG3) 

TWRDl  =  INTGRL  (TWRDl  ,DRT1  ,DELT) 
TWRD2  =  INTGRL  (TWRD2  ,DRT2  ,DELT) 
TWRD3  =  INTGRL  (TWRD3  ,DRT3  ,DELT) 
TWRDl  =  AMAXl  (0.0,  TWRDl) 

TWRD2  =  AMAXl  (0.0,  TWRD2) 

TWRD3  =  AMAXl  (0.0,  TWRD3) 


♦  - ^IfREMOBl  equals  zero  and  TRANSl  equals  zero  and  DVS  greater  that  one,  all  biomass  of 

*  cohort  1  is  added  to  Cohort  2.  Therefore  all  biomass  in  cohort  1  is  gone  &  shouldn't  come  back 
IF  (DVS.GT.1.0  .AND.  REMOBl.EQ.0.0  .AND.  TRANS1.EQ.0.0)THEN 

TWLG2  =  TWLG2  +  TWLGl 
TWRG2  =  TWRG2  +  TWRGl 
TWSG2  =  TWSG2  +  TWSGl 
TWLGl  =  0.0 
TWSGl  =  0.0 
TWRGl  =  0.0 
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ENDIF 


♦  - ^If  REMOB2  equals  zero  and  TRANS2  equals  zero  and  DVS  greater  that  two,  all  biomass  of  cohort 

*  2  is  added  to  cohort  3.  Therefore  all  biomass  in  cohort  2  is  gone  &  shouldn't  come  back 
IF(DVS.GT.2.0.AND.REMOB2.EQ.0.0.AND.TRANS2.EQ.0.0.AND.LAT.LE.33) 

&THEN 

TWLG3  =  TWLG3  +  TWLG2 
TWRG3  =  TWRG3  +  TWRG2 
TWSG3  =  TWSG3  +  TWSG2 
TWLG2  =  0.0 
TWSG2  =  0.0  . 

TWRG2  =  0.0 
ENDIF 

♦  - ^Total  plant  weights 

TWLVG  =  TWLGl  +  TWLG2  +  TWLG3 
TWLVD  =  TWLD 1  +  TWLD2  +  TWLD3 
TWRTG  =  TWRGl  +  TWRG2  +  TWRG3 
TWRTD  =  TWRDl  +  TWRD2  +  TWRD3 
TWSTG  =  TWSGl  +  TWSG2  +  TWSG3 
TWSTD  =  TWSD 1  +  TWSD2  +  TWSD3 

ELSE  IF  (ITASK.EQ.4)  THEN 

*  TERMINAL  SECTION 


♦  - ^Terminal  calculations 

♦  - ^Terminal  output 

CLOSE  (lUNITD) 
ENDIF 

ITOLD  =  ITASK 

RETURN 

END 
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♦ _ 


3,1  ASTRO 


♦** 


SUBROUTINE  ASTRO 
Authors:  Daniel  van  Kraalingen 
Date  :  9  August  1987 
Modified  by  Jan  Goudriaan  4  Febr  1988 
*  Modified  by  Jan  Goudriaan  and  Kees  Spitters  7  December  1989 

Purpose:  This  subroutine  calculates  astronomic  day  length  and  photoperiodic  day  length 
(see  CABO-TPE  report  #?)  and  diurnal  radiation  characteristics  such  as  daily  integral  of 
sine  of  solar  elevation,  solar  constant.  Measured  daily  total  of  global  radiation  is  used  to 
find  atmospheric  transmissivity  and  fraction  diffuse  radiation 
FORMAL  PARAMETERS:  (I=input,0=output,C=control,IN=init,T=time) 


name  meaning 

DAY  Day  number  (Jan  1st  =  1) 

LAT  Latitude  of  tlie  site 

DTR  Measured  daily  total  global  radiation 

SC  Solar  constant 

DSO  Daily  extraterrestrial  radiation 

SINLD  Seasonal  offset  of  sine  of  solar  height 

COSLD  Amplitude  of  sine  of  solar  height 

DAYL  Astronomical  day  lengtli  (base  =  0  degrees) 

DSINB  Daily  total  of  sine  of  solar  height 

DSINBE  Daily  total  of  effective  solar  height 


units 


degrees 
Jm-2d-l 
J  m-2  s-1 
Jm-2d-l 


h 

s 

s 


*  FATAL  ERROR  CHECKS  (execution  terminated,  message) 

*  condition 

*  LAT  >67,  LAT  <-67 

* 

*  SUBROUTINES  and  FUNCTIONS  called  :  none 

* 

*  FILE  usage  :  none 

♦  _ — _ — _ _ _ 


class 


♦ 

♦ 

—  * 


I 

I 

I 

O 

O 

o 

o 

o 

o 

o 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

— * 


SUBROUTINE  ASTRO  (DAY,LAT,SC,DSO,SINLD, COSLD, 
$  DAYL,DSINB, DSINBE) 

IMPLICIT  REAL  (A-Z) 


* - ^PI  and  conversion  factor  from  degrees  to  radians 

PARAMETER  (PI=3. 141592654,  RAD=0.017453292) 


* - Check  on  input  range  of  parameters 

IF  (LAT.GT.67.)  STOP  BRROR  IN  ASTRO:  LAT  >  67 
IF  (LAT.LT.-67.)  STOP  ’ERROR  IN  ASTRO:  LAT  <-67’ 


* - ^Declination  of  the  sun  as  fiinction  of  daynumber  (DAY) 

DEC  =  -ASIN(SIN(23.45*RAD)*COS(2.*PP(DAY+10.)/365.)) 


* - SINLD,  COSLD  and  AOB  are  intermediate  variables 

SINLD  =  SIN(RAD*LAT)*SIN(DEC) 

COSLD  =  COS(RAD*LAT)*COS(DEC) 

AOB  =SINLD/COSLD 
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*  - ^Daylength  (DAYL) 

DAYL  =  12.0*(1.+2.*ASIN(AOB)/PI) 

DSINB  =  3600.*(DAYL*SINLD+24.*COSLD*SQRT(1.-AOB*AOB)/PI) 
DSINBE=3600.*(DAYL*(SINLD+0.4*(SINLD*SINLD+COSLD*COSLD*0.5))+ 
$  12.0*COSLD*(2.0+3.0*0.4*SINLD)*SQRT(l.-AOB*AOB)/PI) 

*  - Solar  constant  (SC)  and  daily  extraterrestrial  (DSO) 

SC  =  1370.*(1.+0.033*COS(2.*PI*D  AY/365.)) 

DSO  =  SC*DSINB 

RETURN 

END 
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♦**  3,2TOTASS  *** 

♦ _ * 

*  SUBROUTINE  TOTASS  * 


*  Authors:  Daniel  van  Kraalingen  * 

*  Date  :  10  December  1987  * 

*  Modified  by  Jan  Goudriaan  5  February  1988  * 

*  Modified  by  Jan  Goudriaan  and  Kees  Spitters  7  December  1989  * 

*  Units  modified  by  Elly  Best  &  Will  Boyd  28  July  1995  * 

*  Purpose:  This  subroutine  calculates  daily  total  gross  assimilation  (DTGA)  by  performing  a  * 

*  Gaussian  integration  over  time.  At  three  different  times  of  the  day,  radiation  is  computed  * 

*  and  used  to  determine  assimilation  whereafter  integration  takes  place.  (Source:  Post-graduate  * 

*  Course  ‘Simulation  of  plant  growth  and  crop  production.  Pontignano,  Siena,  Italy;  3-12  * 

*  November,  1992.  Dept.  Theor.  Production  Ecol.  (TPE-WAU),  Wageningen  Agricultural  * 


*  University,  and  DLO-Centre  for  Agrobiological  Research  (CABO-DLO).)  * 

♦  ♦ 

*  FORMAL  PARAMETERS:  (I=input,0=output,C=control,IN=init,T=time)  * 

*  name  meaning  units  class  * 

*  SC  Solar  constant  J  ni-2  s-1  I  * 

*  D AYL  Day  length  (base  =  0  degrees)  h  I  * 

*  SINLD  Intermediate  variable  in  calculating  solar  declination  -  I  * 

*  COSLD  Intermediate  value  in  calculating  solar  height  -  I  * 

*  DSINBE  Daily  total  of  effective  solar  height  s  I  * 

*  DTR  Measured  daily  total  of  global  radiation  J/m2/d  I  * 

*  RC  Reflection  coefficient  of  irradiation  at  water  surface  (relative)-  I  * 

*  L  Water  type  specific  light  extinction  coefficient  -  I  * 

*  K  Plant  species  specific  light  extinction  coefficient  -  I  * 

*  AMAX  Assimilation  rate  at  light  saturation  for  individual  shoots  gC02/gDW/h  I  * 

*  EE  Initial  light  use  efficiency  for  individual  shoots  g  C02/J  I  * 

*  TL  Thickness  per  plant  layer  m  I  * 

*  DEPTH  Water  depth  ml* 

*  RCHSHST  Relation  coefficient  shoot  weight-stem  length  m/g  DW  I  * 

*  TGW  Total  live  plant  dry  weight  g  DW/ni2  I  * 

*  FGROS  Instantaneous  assimilation  rate  of  whole  canopy  g  C02/  m2  soil/h  O  * 

*  FL  Leaf  diy  matter  allocation  to  each  layer  of  plant  -  I  * 

*  WLV  Diy  weight  of  leaves  gDW/m2  I  * 

*  WST  Diy  weight  of  stems  gDW/m2  I  * 

*  HAR  Harvesting  -  I  * 

*  HARD  AY  Harvesting  day  number  d  I  * 

*  HARDEP  Harvesting  deptli  m  I  * 

*  DTGA  Daily  total  gross  assimilation  gC02/m2/d  O  * 

«  * 

*  SUBROUTINES  and  FUNCTIONS  called  :  ASSIM  * 

*  * 

*  FILE  usage  :  none  * 

♦  _ * 


SUBROUTINE  TOTASS  (SC,DAYL,SINLD,COSLD,DSINBE,DTR,RC,L,K, 
$  AMAX,EE,TL,DEPTH,RCSHST,TGW,FGROS,FL, 

$  WLV,WST,DAY,HAR,HARDAY,HARDEP,DTGA,IRS) 

IMPLICIT  REAL(A-Z) 

REAL  XGAUSS(3),  WGAUSS(3) 
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INTEGER  II,  IGAUSS 
PARAMETER  (PI=3. 14 1592654) 

DATA  IGAUSS  73/ 

DATA  XGAUSS  70. 1 127,  0.5000,  0.88737 
DATA  WGAUSS  70.2778,  0.4444,  0.27787 


*  - Assimilation  set  to  zero  &  three  different  times  of  the  day  (HOUR) 

DTGA  =  0. 

DO  10  U=1,IGAUSS 

♦  - At  the  specified  HOUR,  radiation  is  computed  and  used  to  compute  assimilation 

HOUR=  12.0+DAYL*0.5*XGAUSSai) 


* - Sine  of  solar  elevation 

SINE  =  AMAX1(0.,SINLD+COSLD*COS(2.*PI*(HOUR+12.)724.)) 


♦  - ^Diffuse  light  fraction  (FRDIF)  from  atmospheric  transmission  (ATMTR) 

P7^  =0.5*DTR*SINB*(1.+0.4*SINB)7DSINBE 
ATMTR  =  P7U17(0.5*SC*SINB) 

FRDBF  =  1.47-1.66*ATMTR 

IF  (ATMTR.LE.0.35.AND.ATMTR.GT.0.22)  FRDIF=l.-6.4*(ATMTR-0.22)**2 
IF  (ATMTR.LE.0.22)  FRDIF=1. 

FRDIF  =  AMAX1(FRDIF,0.15+0.85*(1.-EXP(-0.17SINB))) 

♦  - Diffuse  PAR  (PARDIF)  and  direct  PAR  (PARDIR) 

PAR  =  0.5*DTR*SINB*(1.+0.4*SINB)7DSINBE 
PARDIF  =  MIN(PAR,SINB*FRDIF*ATMTR*0.5*SC) 

PARDIR  =  PAR-PARDIF 

CALL  ASSIM 

$(PARDIR,PARDIF,RC,L,K,AMAX,EE,TL,DEPTH,RCSHST,TGW, 

$FL,WLV,WST,DAY,HAR,HARDAY,HARDEP,II,FGROS,IRS) 

♦  - ^Integration  of  assimilation  rate  to  a  daily  total  (DTGA) 

DTGA  =  DTGA+FGROS*WGAUSS(II) 

10  CONTINUE 

DTGA  =  DTGA*DAYL 

RETURN 

END 
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**♦  3.3  ASSIM  *** 

♦ _ _ * 

♦  * 

*  Authore:  Elly  Best*  Will  Boyd  * 

*  Date  :  16  July  1997  * 

*  Modified  by  Jan  Goudriaan  5  Februaiy  1988  * 


*  Purpose:  This  subroutine  performs  an  instantaneous  calculation  of  light  profile  in  the  water  * 

*  column,  light  absorbed  by  the  plant  tissue  available  for  photosynthesis,  and  assimilation  at  * 

*  all  tliese  depth  layers.  The  depth-integrated  variable  is  FGROS.  At  harvesting,  the  plant  * 


*  material  is  removed  per  depth  layer  from  the  existing  biomass. 

* 

* 

* 

*  FORMAL  PARAMETERS:  (I=inpuLO=output,C=control,IN=init,T=time) 

* 

*  name 

*  — — 

meaning 

units  class  * 

.  _ _  * 

*  PARDIR 

Instantaneous  flux  of  direct  radiation  (PAR) 

W/m2 

I  * 

*  PARDIF 

Instantaneous  flux  of  diffuse  radiation(PAR) 

W/m2 

I  * 

*  RC 

Reflection  coefficient  of  irradiation  at  water  surface  (relative) 

- 

I  * 

*  L 

Water  type  specific  light  extinction  coefficient 

m-I 

I  * 

*  K 

Plant  species  specific  light  extinction  coefficient 

m2/gDW 

I  * 

*  AMAX 

Assimilation  rate  at  light  saturation  for  individual  shoots 

gC02/gDW/h 

I  * 

*  EE 

Initial  light  use  efficiency  for  individual  shoots 

g  C02/J 

I  * 

*  TL 

Thickness  per  plant  layer 

m 

I  * 

*  DEPTH 

Water  dep^ 

m 

I  * 

*  RCHSHST 

Relation  coefficient  shoot  weight-stem  lengtli 

m/gDW 

I  * 

*  TGW 

Total  live  plant  dry  weight 

gDW/m2 

I  * 

*  FL 

Leaf  dry  matter  allocation  to  each  layer  of  plant 

- 

I  * 

*  WLV 

Dry  weight  of  leaves 

gDW/m2 

I  * 

*  WST 

Dry  weight  of  stems 

gDW/m2 

I  * 

*  HAR 

Harvesting 

- 

I  * 

*  HARDAY 

Harvesting  day  number 

d 

I  * 

*  HARDER 

Harvesting  depth 

m 

I  * 

*  II 

Counter  in  DO  LOOP,  indicates  I  of  3  times  per  day  (HOUR) 

- 

I  * 

*  FGROS 

Instantaneous  assimilation  rate  of  the  crop 

g  C02/m2/h 

O  * 
* 

*  SUBROUTINES  called  :  none 

* 

*  FUNCTIONS  called :  AFGEN 

3k 

* 

* 

*  FILE  usage : 

:  none 

* 

* 

SUBROUTINE  ASSIM  (PARDIR,PARDIF,RC,L,K,AMAX,EE,TL, 
$  DEPTH,RCSHST,TGW,FL,WLV,WST,DAY, 

$  HAR,HARDAY,HARDEP,II,FGR0S,IRS) 

IMPLICIT  RE  AL(A-Z) 

REAL  DMPC(6),  SC(IOO),  IRZ(100) ,  IABS(1(K)),  lABSL(lOO) 
REAL  fflG(lOO),  AH(I00),  REDF(IOO),  SumZ,  BotBio 
INTEGER  IMNl,  IRED,  I,  LOOP,  Layers,  LBelow,  ILAY,  II 
PARAMETER  (IMNI  =  40) 

REAL  REDFT(IMNI),  DMPCT(IMNI) 


* - ^Read  AFGEN  functions 

CALLED  AREA  (■REDFT’,REDFT, IMNI  ,IRED  ) 
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CALL  RDAREA  CDMPCT  ',DMPCT,  IMNl  ,ILAY  ) 

*  - Irradiation  just  beneath  the  water  surface 

IRS  =  PARDIR  +  PARDIF 
IRZ(1)  =  IRS*(1.0-RC) 

*  - Canopy  assimilation  is  set  to  zero 

FGROS  =  0. 

*  - Calculate  stem  lengtli 

STEMLE  =  AMINl(Depth+.0995,  {RCSHST*(WLV+WST))) 

IF  (STEMLE  .GT.  Depth+.08)THEN 

*  - Determine  total  number  of  layers  in  the  given  water  depth 

LOOP  =  INT  (Depth/TL  +  0. 1)  +  1 

*  - ^LOOP  should  never  be  less  tlan  6  since  DEPTH  shouldn't  be  less  than  .5m 

IF  (LOOP  .LE.  5)  LOOP  =  6 

*  - Distribute  61%  of  total  plant  biomass  in  1st  5  layers 

DO  10 1  =  1,5 
VAL  =  REAL(I) 

DMPC(I)  =  LINT  (DMPCT,ILAY,VAL) 

SC(I)  =  TGW  ♦  DMPC(I) 

10  CONTINUE 

*  - If  water  depth  is  at  least  Im  -  use  METHOD  1  for  distribution  of 

*  biomass  beyond  1st  5  layers;  otherwise,  use  METHOD2 

*  METHOD  1 

If  (LOOP  .GE.  10)THEN 

*  - Distribute  39%  of  biomass  in  tlie  lower  layers  (including  last  layer) 

*  with  biomass  gradually  decreasing  toward  tlie  bottom 

*  LOOP  (integer) ..  Number  of  0.  Im  water  layers 

*  LAYERS  (integer)  ..  Layers  remaining  after  initial  5 

*  - SUMZ  (real) ..  Summation  of  layers  6  tluough  LOOP 

*  - ^LBELOW  (integer) ..  Layer  number  going  from  bottom  to  top 

*  - 5  in  next  statement  is  for  tite  1st  5  layers  of  the  plant 

LAYERS  =  LOOP  -  5 

SUMZ  =  (LAYERS/2.0)  *  (LAYERS+1.0) 

DO  20 1  =  6,LOOP 

LBELO W  =  LAYERS  -  (1-5)  +1.0 

SC(I)  =  (LBELOW/SUMZ)  *  (TGW  *  0.39) 

20  CONTINUE 

*  - METHOD  2 

ELSE 

*  - If  water  depth  is  less  than  Im  -  put  5%  of  total  biomass  in  each  layer 

*  remaining  -  subtract  fi’om  the  39%  biomass  reserved  for  lower  layers 
BotBio  =  TGW  *  0.39 
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DO  21 1  =  6, LOOP 
SC(I)  =  TGW  *  0.05 
BotBio  =  BotBio  -  SC(I) 

21  CONTINUE 

*  - ^Redistribute  difference  "BotBio"  over  the  top  5  layers  proportionally 

DO  22 1  =  1,5 

sea)  =  sea)  +  (DMPCa)*BotBio) 

22  CONTINUE 

ENDEF 

*  - ^Harvesting 

IF  OiAR  .EQ.  1.  .AND.  DAY  .EQ.  HARDAY)THEN 
IF  OlARDEP  .GT.  DEPTH)  HARDEP  =  DEPTH 
D0  25I=  l,HARDEP/.l 
SCO)  =  0.0 
25  CONTINUE 

*  - ^Reset  total  live  weight  (TGW)  to  zero 

IFai.EQ.  1)TGW  =  0.0 
ENDIF 

DO  50 1  =  l,LOOP 


♦  - ^Total  irradiation  on  top  of  stratum  I 

mza+l)  =  IRZ(I)  *  EXP(-TL*  L  -  K*  SC(I)) 

IF(SC(I)  .EQ.  0.0)  GOTO  30 

*  - ^Radiation  absorbed  by  macrophyte  community 

lABSa)  =  (IRZ(I)-IRZ(I+l))*SC(I)*K/aC*SC(I)  +  TL*L) 


*  - ^Radiation  absorbed  by  leaves,  excluding  bottom  layer 

IFa  .LT.  LOOP)  IABSL(I)  =  lABS(I)  *  FL 
IFaABSL(I)  .EQ.  0.0)GOTO  30 

*  - ^Height  on  top  of  stratum  I  measured  from  the  water  surface 

HIGa)  =  TL  *  O-OOP  - 1) 

*  - Absolute  height  of  vegetation  on  top  of  stratum  I,  measured 

*  from  the  top  of  the  plant 
AHa)  =  STEMLE  -  HIGa) 

*  - ^Reduction  factor  over  the  vertical  of  tlie  vegetation 

REDF(I)  =  LINT(REDFT,IRED,AHa)) 

*  - Instantaneous  C02  assimilation  rate  per  deptli  layer 

FGL  =  SCa)*AMAX*REDF(I)*(l.-EXP(-EE*IABSLa)*3600.  / 
$  (AMAX*REDF(I)*SC(I)))) 

GOTO  40 
30  FGL  =  0.0 

40  FGROS  =  FGROS  +  FGL 

*  - If  plants  are  liarvested,  live  plant  weight  is  recalculated 
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IF  (HAR.EQ.1  .AND.  DAY.EQ.HARDAY  .AND.  lI.EQ.l)  THEN 
TGW  =  TGW  +  SC® 

ENDBF 

50  CONTINUE 
ENDIF 

RETURN 

END 

^![t****************************************************************************** 

***  3.4  CHORT2  ♦** 


*  Authors:  Elly  Best  &  Will  Boyd 

*  Date  :  16  July  1997 

*  Purpose:  This  subroutine  describes  the  behavior  of  the  second  plant  cohort 

*  _ 


♦ 

* 

* 

* 

.* 


SUBROUTINE  CHORT2(AMAX  ,AMTMP  ,AMTMPT,AMX  ,ASRQ  ,COSLD,CRRIZ, 
&  CVT  ,DAVTMP,DAY  ,DAYL  ,DDTMP  ,DELT  .DEPTH , 

&  DLV2  ,DRT2  .DSINB  ,DSINBE,DST2  ,DTEFF,DTGA  , 

&  DVR  ,DVS  ,EE  ,FGROS,FL  ,FLT  ,FLV  , 

&  FLVT  ,FRT  ,FRTT  ,FST  ,FSTT  ,GLV  .GPHOT , 

&  GRT  ,GST  ,GTW  ,HAR  ,HARDAY,HARDEP,IKT  , 

&  ILAMTM,ILFLT  ,ILFLVT,ILFRTT,ILFSTT,ILRDRT,ILRDST, 

&  ILT  ,K  ,KT  ,L  ,LAT  ,LT  , 

&  MAINRT,MAINT  ,NGLV2  ,NGRT2  ,NGST2  ,RC  , 

&  RCSHST,RDD  ,RDR  .RDRIZ  ,RDRT  ,RDST  .REDAM , 

&  REMOBl.REMOB2,ROC  ,SC  .SINLD  .TBASE  ,TEFF  , 

&  TGRIZ.TGW  .TIME  .TL  .TRAFAC.TRANS1.TRANS2. 

&  TREMOB.TW  .TWLDl  .TWLD2  .TWLGl  .TWLG2  .TWRDl . 

&  TWRD2  .TWRG 1 ,  TWRG2  .TWGRIZ.TWRIZD.TWSD  1 . 

&  TWSD2  .TWSGl.  TWSG2  .WLV  .WST  .IRS) 

IMPLICIT  REAL  (A-Z) 


*  - ^Formal  parameters 

INTEGER  IKT  .ILAMTM.ILFLT  .ILFLVT.ILFRTT.ILFSTT 
INTEGER  ILRDRT.ILRDST.ILT 
LOGICAL  OUTPUT,  TERMNL 

IF  ®VS  .LT.  1.0)  GOTO  100 

* — Call  to  SBRT  ASTRO  to  introduce  day  length  into  MAIN  for  tentative 

*  relationship  REMOB2-DAYL;  otherwise  this  call  can  be  made  at  ♦ASTROOO 
CALL  ASTRO 

$®AY,LAT,SC,DSO,SINLD,COSLD,DAYL,DSINB,DSINBE) 


* - Calculation  of  diy  matter  and  its  partitioning  over  the  plant  organs 

TW  =  TGW  +  (TWLD1+TWSD1+TWRD1+TWLD2+TWSD2+TWRD2) 
FLV  =  LINT{FLVT  ,ILFLVT,DVS) 

FST  =  LINT(FSTT  .ILFSTT.DVS) 

FRT  =  LINT(FRTT  ,ILFRTT,DVS) 

FL  =LINT(FLT  ,ILFLT,DVS) 


Appendix  A  Model  Listing 


A21 


*  - Growth  of  plant  organs,  maintenance  respiration  and  translocation 

*  Calculation  assimilate  requirement  for  plant  organ  formation. 

ASRQ  =  1.54 

MAINTS  =  0.016*TWLG2+0.01*TWSG2+0.015*TWRG2 
MAIN!  =MAINTS*TEFF 
MAINRT  =  0.005  *  TWGRIZ  *  TEFF 

*  - Carbohydrate  behavior:  remobilization  from  rliizomes  for  plant  formation  at  proper  day  length 

*  and  temperature  conditions  (presently  related  to  DVS);  translocation  from  plants  to  form 

*  rhizomes,  provided  plants  are  present 
TRANS2  =  0.0 

TREMOB  =  0.0 

TWGRIZ  =  AMAXl  (CRRIZ,  TWGRIZ) 

TGRIZ  =  TWGRIZ 

TWRIZD  =  INTGRL  (TWRIZD,  RDRIZ  ,  BELT) 

IF  (DVS.GE.1.0  .AND.  DVS.LT.1.63)  THEN 
IF  (GPHOT  .LT.  MAINT)  THEN 
IF  (TWGRIZ  .GT.  CRRIZ)  THEN 
TREMOB  =  INTGRL  (TREMOB,  REMOB2,  BELT) 

REMOB2  =  ROC  *  TWGRIZ 
ELSE 

WRITE(*,*)'  Vegetation  is  dying  * 

REMOB2  =  0.000001 
ENDIF 

ENDIF 

ELSE 

REMOB2  -  0.000001 
ENDIF 

*  - If  the  carbohydrates  fill  the  rhizomes  by  TRANS  1  from  cohort  1,  REMOB  1  stays  zero;  TRANS2 

*  remains  zero,  if  tlie  carbohydrates  leave  the  rhizomes  by  REMOB2  to  fonn  the  plants  of  cohort 

*  but  also  if  no  plant  biomass  of  cohort2  is  present.  The  plants  of  cohort2  are  formed  by  REMOB2 

*  only  at  a  certain  DVS 

IF  (TRANS  1  .GT.  0.0)REMOB1  =  0.0 
IF  (DVS  .GE.  2.0)REMOB1  -  0.0 
IF  (DVS  .GE.  2.0)TRANS1  =  0.0 
IF  (REMOB2  .GT.  0.00000 1)TRANS2  =  0.0 

*  - ^Relative  death  rates 

RDR  =  INSW  (DVS-2.001,0.,  LINT(RDRT,ILRDRT,DDTMP)) 

RDS  =  INSW  (DVS-2.001,0.,  LINT(RDST,ILRDST,DDTMP)) 


♦ - Development  rates 

IF(DAVTMP  .LT.  3.0)  THEN 
DVR  =  0.0 

ELSE  IF  (DVS.LE.l.)  THEN 
DVR  =  0.022*DAVTMP/30 
ELSE  IF  PVS.LE.6.0)  THEN 
DVR  =  0.015*DAVTMP/30 
ENDIF 
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•Calculation  of  effective  daytime  temperature 
DTEFF  =  AMAX1(0.,DAVTMP-TBASE) 


•Calculation  of  dead  plant  material  per  organ 
DLV2  =  TWLG2  *  RDR 
DST2  =  TWSG2  *  RDR 
DRT2  =  TWRG2  *  RDR 


•Shoot  photosynthesis  at  light  saturation  and  daytime  temperature  effect  on  shoot  photosynthesis 
AMAX  =  AMAX1(0.00001,AMX  *  AMTMP) 

AMAX  =  AMAX*  REDAM 

AMTMP  =  LINT(AMTMPT,a.AMTM,DDTMP) 


*  - ^Before  calling  TOTASS,  detennine  light  extinction  coefficients  of  plants  (K)  and  of  water  (L) 

L  =  LINT(LT,ILT,TIME) 

K  =  LINT(KT,IKT,DVS) 

*  - ^Daily  total  gross  assimilation 

CALL  TOTASS 

$  (SC,DAYL,SINLD,COSLD,DSINBE,RDD,RC,L,K,AMAX,EE, 

$  TL,DEPTH,RCSHST,TGW,FGROS,FL,WLV,WST, 

$  DAY,HAR,HARDAY,HARDEP,DTGA,IRS) 


*  - If  liarvesting  takes  place,  weights  various  plant  organs  must  be  recalculated 

*  (TWLVG,TWSTG,TWRTG,TW) 

IF(HAR  .EQ.  1.  AND.  DAY  .EQ.  HARD  AY)  THEN 
TWLG2  =  FLV  *  TGW 
TWSG2  =  FST  *  TGW 
TWRG2  =  FRT  *  TGW 

TW  =  TGW  +  (TWLD1+TWSDH-TWRD1+TWLD2+TWSD2+TWRD2) 
ENDIF 


* - Conversion  assimilated  C02  to  CH20 

GPHOT  =  DTGA  *  30./44. 


*  - Induction  of  flowering  at  DVS=2;  flowering  occurs  10  days  after  induction. 

*  Induction  of  flowering,  translocation  and  senescence  occur  simultaneously  . 
IF  (DVS  .GE.  I.63)THEN 

IF  ((TWLG2+TWSG2+TWRG2)  .GT.  0.0)THEN 

IF  (TGW  .GE.  TWGRIZ)  THEN 
IF  (TRANS  1  .EQ.  0.0)  THEN 
TRANS2  =  CVT  *  GPHOT*((TWLG2+TWSG2+TWRG2)/ 

$  (TWLG1+TWSG1+TWRGI+TWLG2+TWSG2+TWRG2))  *  TRAFAC 

ENDIF 

ENDIF 

ENDIF 

*  - If  there  is  no  plant  biomass  REMOB2  must  stay  at  zero 

IF  ((TWLG2+TWSG2+TWRG2)  .EQ.  0.0)REMOB2  =  0.0 

ENDIF 
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*  - When  TRANS2  gets  a  value  set  REM0B2  to  0 

If  (TRANS2  .GT.  0.0)REMOB2  =  0.0 

*  - The  total  weight  of  the  green  rhizomes  increases  by  translocation  (TRANS2  from  CHORT2) 

*  and  decreases  by  remobilization  and  rhizome  maintenance  respiration 
TWGRIZ  =  INTGRL  (TWGRIZ,  -(((-TRANS2  +REMOB2+  MAINRT)/1.242) 

$  +(RDRIZ*TGRIZ)),  DELT) 

*  - Total  and  net  growth  rates 

GTW  =  ((REMOB2*CVT)  +  GPHOT  -  TRANS2  -  MAINT)  /  ASRQ 
GRT  =  FRT  *  GTW 
GST  =  FST  *  GTW 
GLV  =  FLV*GTW 

NGLV2  =  GLV  -  DLV2 
NGST2  =  GST  -  DST2 
NGRT2  =  GRT-DRT2 


100  RETURN 
END 

***  3.5  CHORT3  *** 


* 

*  Authors:  Elly  Best  &  Will  Boyd 

*  Date  :  18  August  1997 

*  Purpose:  This  subroutine  describes  the  behavior  of  the  third  plant  cohort 

*  _ _ _ 


* 

* 

* 

* 

.* 


SUBROUTINE  CHORT3(AMAX  ,AMTMP  ,AMTMPT,AMX  ,ASRQ  ,COSLD,CRRIZ, 
&  CVT  ,DAVTMP,DAY  ,DAYL  ,DDTMP  ,DELT  .DEPTH  , 

&  DLV3  ,DRT3  ,DSINB  ,DS1NBE,DST3  .DTEFF  ,DTGA  , 

&  DVR  ,DVS  ,EE  ,FGROS,FL  ,FLT  ,FLV  , 

&  FLVT  ,FRT  ,FRTT  ,FST  ,FSTT  ,GLV  , GPHOT  , 

&  GRT  ,GST  ,GTW  ,HAR  ,HARDAY,HARDEP,IKT  , 

&  ILAMTM,ILFLT  ,1LFLVT,1LFRTT,1LFSTT,1LRDRT,1LRDST, 

&  ILT  ,K  ,KT  ,L  ,LAT  ,LT  , 

&  MAINRT,MA1NT  ,NGLV3  ,NGRT3  ,NGST3  ,RC  , 

&  RCSHST,RDD  ,RDR  ,RDRIZ  ,RDRT  ,RDST  , RED  AM , 

&  REMOBl,REMOB3,ROC  ,SC  ,S1NLD  ,TBASE  ,TEFF  , 

&  TGRIZ.TGW  .TIME  .TL  ,TRAFAC.TRANS1.TRANS2. 

&  TRANS3.TREMOB.TW  .TWLDl  .TWLD2  ,TWLD3  .TWLGl . 

&  TWLG2  .TWLG3  .TWRDl  .TWRD2  .TWRD3  .TWRGl  .TWRG2  . 

&  TWRG3  .TWGRIZ.TWRIZD.TWSD1  .TWSD2  .TWSD3  . 

&  TWSGl.  TWSG2  .TWSG3  .WLV  .WST  ,IRS) 

IMPLICIT  REAL  (A-Z) 


* - Formal  parameters 

INTEGER  IKT  .ILAMTM.ILFLT  .ILFLVT.ILFRTT.ILFSTT 
INTEGER  ILRDRT.ILRDST.ILT 
LOGICAL  OUTPUT.  TERMNL 
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IF  (DVS  XT.  2.0)  GOTO  100 

♦-—Call  to  SBRT  ASTRO  to  introduce  day  length  into  MAIN  for  tentative 

*  relationship  REMOB3-DAYL;  otherwise  this  call  can  be  made  at  *ASTROOO 

CALL  ASTRO 

$  (dayxat,sc,dso,sinld,cosld,dayl,dsinb,dsinbe) 

*  - Calculation  of  dry  matter  and  its  partitioning  over  the  plant  organs 

TW  =  TGW+(TWLD1+TWSD1+TWRD1+TWLD2+TWSD2+TWRD2+TWLD3+TWSD3+TWRD3) 
FLV  =  LINT(FLVT  ,ILFLVT,DVS) 

FST  =  LINT(FSTT  ,ILFSTr,DVS) 

FRT  =  LINT(FRTT  ,ILFRTT,DVS) 

FL  =LINT(FLT  ,ILFLT,DVS) 


*  - Growth  of  plant  organs,  maintenance  respiration  and  translocation 

*  Calculation  assimilate  requirement  for  plant  organ  formation. 

ASRQ  =1.54 

MAINTS  =  0.016*TWLG3+0.01*TWSG3+0.015*TWRG3 
MAINT  =  MAINTS*TEFF 
MAINRT  =  0.005  *  TWGRIZ  *  TEFF 

*  - Carbohydrate  behavior;  remobilization  from  rhizomes  for  plant  formation  at  proper  day  length  and 

*  temperature  conditions  (presently  related  to  DVS);  translocation  from  plants  to  form  rhizomes, 

*  provided  plants  are  present 
TRANS3  =  0.0 
TREMOB  =  0.0 

TWGRIZ  =  AMAXl  (CRRIZ,  TWGRIZ) 

TGRIZ  =  TWGRIZ 

TWRIZD  =  INTGRL  (TWRIZD,  RDRIZ  ,  DELT) 

IF  (DVS.GE.2.0  .AND.  DVS.LT.2.447)  THEN 
IF  (GPHOT  .LT.  MAINT)  THEN 
IF  (TWGRIZ  .GT.  CRRIZ)  THEN 
TREMOB  =  INTGRL  (TREMOB,  REMOB3,  DELT) 

REMOB3  =  ROC*  TWGRIZ 
ELSE 

WRITE(*,*)'  Vegetation  is  dying ' 

REMOB3  =  0.000001 
ENDIF 

ENDIF 

ELSE 

REMOB3  =  0.000001 
ENDIF 

*  - If  tlie  carbohydrates  fill  the  rhizomes  by  TRANS  1  from  cohort  1,  REMOB  1  stays  zero;  TRANS2 

*  remains  zero,  if  the  carbohydrates  leave  the  rhizomes  by  REMOB3  to  form  the  plants  of  cohort3, 

*  but  also  if  no  plant  biomass  of  cohort2  is  present.  The  plants  of  cohort2  are  formed  by  REMOB3 

*  only  at  a  certain  DVS 

IF  (TRANSl  .GT.  0.0)REMOB1  =  0.0 
IF  (DVS  .GE.  2.0)REMOB1  =  0.0 
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IF  (DVS  .GE.  2.0)TRANS1  =  0,0 
IF  (DVS  .GE.2.3)TRANS2  =  0.0 
IF  (REMOB3  .GT.  0.00000 1)TRANS3  =  0.0 


* - Relative  death  rates 

RDR  ==  INSW  (DVS-3.5001,0.,  LINT(RDRT,ILRDRT,DDTMP)) 
RDS  =  INSW  (DVS-3.5001,0.,  LINT(RDST,ILRDST,DDTMP)) 


♦  - Development  rates 

IF(DAVTMP  .LT.  3.0)  THEN 
DVR  =  0,0 

ELSE  IF  (DVS.LE.  1.)  THEN 
DVR  =  0.022*DAVTMP/30 
ELSE  IF  (DVS.LE.6.0)  THEN 
DVR  =  0.015*DAVTMP/30 
ENDIF 

♦  - ASTROOO 

♦  - Calculation  of  effective  daytime  temperature 

DTEFF  =  AMAX1(0.,DAVTMP-TBASE) 


•Calculation  of  dead  plant  material  per  organ 
DLV3  =TWLG3*RDR 
DST3  =  TWSG3  *  RDR 
DRT3  =  TWRG3  *  RDR 


♦  - Shoot  photosyntliesis  at  light  saturation  and  daytime  temperature  effect  on  shoot  photosynthesis 

AMAX  =AMAX1(0.00001,AMX*  AMTMP) 

AMAX  =  AMAX  *  REDAM 

AMTMP  =  LINT(AMTMPT,ILAMTM,DDTMP) 

*  - ^Before  calling  TOTASS,  detennine  light  extinction  coefficients  of  plants  (K)  and  of  water  (L) 

L  =  LINT(LT,ILT,TIME) 

K  =  LINT(KT,IKT,DVS) 


* - ^Daily  total  gross  assimilation 

CALL  TOTASS 

$  (sc,dayl,sinld,cosld,dsinbe,rdd,rc,l,k,amax,ee, 
$  tl,depth,rcshst,tgw,fgros,fl,wlv,wst, 

$  DAY,HAR,HARDAY,HARDEP,DTGA,IRS) 


*  - If  harvesting  takes  place,  weights  various  plant  organs  must  be  recalculated 

*  (TWLVG,TWSTG,TWRTG,TW) 

IF(HAR  .EQ.  L  and.  day  .EQ.  hard  AY)  THEN 
TWLG3=FLV*TGW 
TWSG3  =  FST  *  TGW 
TWRG3  =  FRT  *  TGW 

TW  =  TGW+(TWLD1+TWSD14-TWRD1+TWLD2+TWSD2+TWRD2+TWLD3+TWSD3+TWRD3) 
ENDIF 


*  - Conversion  assimilated  C02  to  CH20 

GPHOT  =  DTGA  *  30./44. 

*  - Induction  of  flowering  at  DVS=2.447;  flowering  occurs  10  days  after  induction 
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*  Induction  of  flowering,  translocation  and  senescence  occur  simultaneously 
IF  (DVS  -GE.  2.447)THEN 

IF  ((TWLG3+TWSG3+TWRG3)  .GT.  0.0)THEN 

F  (TGW  .GE.  TWGRIZ)  THEN 
F  (TRANS2  .EQ.  0.0)  THEN 

TRANS3  =  CVT  *  GPHOT  *((TWLG3+TWSG3+TWRG3)/ 
$(TWLG1+TWSG1+TWRG1+TWLG2+TWSG2+TWRG2+TWLG3+TWSG3+TWRG3)) 
$*TRAFAC 

ENDF 

ENDF 

ENDF 

* — -If  there  is  no  plant  biomass  REMOB3  must  stay  at  zero 
F  ((TWLG3+TWSG3+TWRG3)  .EQ.  0.0)REMOB3  =  0.0 

ENDF 

*  - When  TRANS3  gets  a  value  set  REMOB3  to  0 

F(TRANS3  .GT.  0.0)REMOB3  =  0.0 


*  - The  total  weight  of  tlie  green  rhizomes  increases  by  translocation  (TRANS3  from  CHORT3) 

*  and  decreases  by  remobilization  and  rhizome  maintenance  respiration 
TWGRIZ  =  INTGRL  (TWGRIZ,  -(((-TRANS3  +REMOB3+  MAINRT)/1.242) 

$  +(RDRIZ*TGRIZ)),  DELT) 

*  - Total  and  net  growth  rates 

GTW  =  ((REMOB3*CVT)  +  GPHOT  -  TRANS3  -  MAINT)  /  ASRQ 
GRT  =  FRT  *  GTW 
GST  =  FST  *  GTW 
GLV  =  FLV*GTW 
NGLV3  =  GLV-DLV3 
NGST3  =  GST  -  DST3 
NGRT3  =  GRT  -  DRT3 

100  RETURN 
END 
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* 


♦ _ _ _ _ _ _ _ _ _ _ _ _ _ _ * 

*MODEL.DATfile  * 

*  contains:  * 

*  -  Initial  constants  as  far  as  specified  with  INCON  statements,  * 

*  -  Model  parameters,  * 

*  -  AFGEN  functions,  * 

*  -  A  SCALE  array  in  case  of  a  general  translation  * 

*  * 

*  File;  MDELMW70.DAT;  to  be  used  as  input  file  for  MILFO.FOR  * 

*  Contains  data  on  biomass  and  water  temperatures  from  Lake  Wingra,  WI,  1970  * 

*  Date:  10-10-97  * 

*  Time:  10:00:06  * 

*  _ _ _ _ _ * 

*  Initial  constants 

CRRIZ  =  12. 

IWGRIZ  =  50. 

IWRIZD  =  0. 

IWLDl  =0. 

IWLD2  =0. 

IWLD3  =0. 

IWLGl  =23.5 
IWLG2  =0. 

IWLG3  =0. 

IWRDl  =0. 

IWRD2  =0. 

IWRD3  =0. 

IWRGl  =3. 

IWRG2  =0. 

IWRG3  =0. 

IWSDl  =0. 

IWSD2  =0. 

IWSD3  =0. 

IWSGl  =23.5 
IWSG2  =0. 

IWSG3  =0. 

NUL  =0. 

REMOBl  =  0. 


*  Model  parameters 

YRNUM  =  1. 
AMX  =  0.0165 
CVT  =  1.05 
DAYEM  =  1. 
DELAY  =  7. 
DEPTH  =  1.5 
EE  =0.000011 
HAR  =0. 
HARDAY=212. 
HARDEP  =  0.8 
NPL  =11. 

RC  =  0.06 
RCSHST  =  12.0 
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RDRIZ  =0.00042 
REDAM  =0.5 
ROC  =0.0576 
TEASE  =  3. 

TL  =0.1 
TRAFAC  =  0.35 


♦  AFGEN  functions 

♦AMD  VST  = 

*  0.001,1., 

*  3.5,1., 

*  5.0,  1. 
AMTMPT  = 

-30.,  0.00001, 

0.,  0.00001, 

5.,  0.18, 

10.,  0.23, 

15.,  0.40, 

20.,  0.63, 

25.,  0.78, 

30.,  0.95, 

35.,  1.0, 

40.,  0.78, 

45.,  0.38, 

50.,  0.00001 
DMPCT  = 

1.0,  .10, 

2.0,  .16, 

3.0,  .17, 

4.0,  .10, 

5.0,  .08 
♦DVRVT  = 

*  -15.,  0., 

*  0.,0., 

*  30.,  0.022 
♦DVRRT  = 

*  -15.,  0., 

*  0.,  0., 

*  30.,  0.015 
FLT  = 

0.,  0.50, 

2.3,  0.50, 

6.0,  0.50 
FLVT  = 

0.,  0.47, 

2.3, 0.47, 

6.0, 0.47 
FSTT  = 

0.,  0.47, 

2.3, 0.47, 

6.0,  0.47 
FRTT  = 

0.,  0.06, 
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2.3,  0,06, 
6.0,  0.06 
KT  = 

0.,  0.006, 
3.5,  0.006, 
6.0,  0.006 
LT  == 

0.,  1.15, 

101.,  1.6, 

117.,  1.4, 

131.,  1.85, 

156.,  1.85, 

173.,  1.8, 

187.,  1.8, 

215.,  1.8, 

243.,  2.0, 

257.,  1.4, 

271.,  1.6, 

299.,  1.4, 

327.,  1.15, 

365.,  1.15 
RDRT  = 
-30.,  0.042, 

0.,  0.042, 

50.. 0.042 
RDST  = 
-30.,  0.042, 

0.,  0.042, 

50.,  0.042 
REDFT  = 

0.0,  1.0, 

1.0,  1.0, 

6.0,  1.0 
TEFFT  = 
-30.,  0.0001, 
0.,  0.0001, 

5.,  1.0, 

20.. 1.0, 

30.,  2.0, 

40.,  4.0, 

45.,  8.0, 

50.,  12.0 
TGWMT  = 

1.,  50., 

141.,  50., 

162.,  120., 

172.,  106., 

192.,  105., 

202.,  130., 

223.,  160., 

233.,  180., 

254.,  220., 

264.,  150., 

365.,  50.7 


A30 


Appendix  A  Model  Listing 


WTMPT  = 

1.,  3.5, 

60.,  3.5, 

100.,  9.0, 

150.,  22.0, 

190.,  25.0, 

220.,  25.0, 

250.,  19.0, 

300.,  9.0, 

340.,  1.6, 

365.,  1.6 
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*  TIMER  file  contains 

*  -  The  used  DRIVER  and  TRACE  in  case  of  GENERAL  translation 

*  -  The  TIMER  variables  used  in  both  translation  modes 

*  -  Additional  TIMER  variables  in  case  of  GENERAL  translation 

*  -  The  WEATHER  control  variables  if  weather  data  are  used 

*  -  Miscellaneous  FSE  variables  in  case  of  FSE  translation 

* 

*  FUe:  MILFO.FOR 

*  Date:  10-10-97 

*  Time:  10:55:06 

*  TIMER  variables  used  in  GENERAL  and  FSE  translation  modes 


*  T _ _ _ _ _ _ _ _ _ _ * 

S  IT'IME  =1.  !  start  time 

FINTIM  =  365.  !  finish  time 

DELT  =1.  !  time  step  (for  Runge-Kutta  first  guess) 

PRDEL  =  1,  !  output  time  step 

IPFORM  =  4  !  code  for  output  table  format: 


!  4  =  spaces  between  columns 
!  5  =  TAB'S  between  colunms  (spreadsheet  output) 

!  6  =  two  column  output 

!  Tlie  string  array  PRSEL  contains  the  output  variables  for  which 
!  fonnatted  tables  have  to  be  made.  One  or  more  times  there  is  a 
!  series  of  variable  names  terminated  by  the  word  <TABLE>. 

!  The  translator  writes  tlie  variables  in  each  PRINT  statement  to 
PRSEL  =  !  a  separate  table. 

*  'DAVTMP', 

*  'DAYL 

*  'DDTMP 

*  'DTEFF 

*  TITGA 
*’DVS 

*  'FGROS 

*  'GPHOT', 

*'IRS 

""MAINRT, 

*  'MAINT', 

*  'REMOBl', 

*  'REMOB2’, 

*  -REMOBB’, 

■TGW', 

*  TGWM', 

*  'TMPSUM', 

*  'TRANSr, 

*  'TRANS2', 

*  TRANSS', 

*  •TW 
•TWGRIZ', 

*  'TWLDl 

*  'TWLD2 

*  'TWLD3 

^  innxT  rr  ■%  i 


♦ 

* 

* 

* 

* 

* 

♦ 
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*  TWLG2 

*  TWLG3 

*  'TWLVD 

*  TWLVG 

*  TWRDl 

*  TWRD2 

*  'TWRD3 

*  'TWRGl 

*  TWRG2 

*  ’TWRG3  *, 

♦TWRIZ’, 

*'TWRIZD’, 

*  'TWRTD 

*  TWRTG 

*  TWSDl  •, 

*  'TWSD2 

*  TWSD3 

*  'TWSGl 

*  'TWSG2 

*  ’TWSG3 

*  TWSTD  *, 

*  'TWSTG 
*'WTMP’, 

•<TABLE>' 

COPINF  =  'N'  !  Switch  variable  whether  to  copy  tlie  input  files 

!  to  tl»e  output  file  ('N'  =  do  not  copy, 

!  'V  =  copy) 

DELTMP  =  'N'  !  Switch  variable  what  should  be  done  with  the 

!  temporary  output  file  ('N'  =  do  not  delete, 

!  •Y'  =  delete) 

IFLAG  =1101  !  Indicates  where  weather  error  and  warnings 

!  go  (1101  means  errors  and  warnings  to  log 
!  file,  errors  to  screen,  see  FSE  manual) 
"“lOESD  =  1991,182  !  List  of  observation  data  for  which  output  is 
!  required.  The  list  should  consist  of  pairs 
!  <year>,<day>  combination 


*  WEATHER  control  variables 


WTRDIR  =  'C:\SYS\WEATHERV 
CNTR  =  'WIS'  !  Country  code 

ISTN  =  1  !  Station  code 

lYEAR  =  1970  !  Year 
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* _ _ _ * 

*  CONTROL.DAT  file  contains:  * 

*  -  File  names  to  be  used  by  FSE  2.1  * 

*  The  input  files  (except  FILEIR)  may  used  in  reruns  * 

*  Up  to  five  input  data  files  may  be  used  (FILEIl-5)  * 

*  _ _ _ * 


FILEON='RES.DAT' 
FILEOL  =  ’MODEL.LOG' 
FILEIR  =  'RERUNS.DAT' 
FILErr  =  ’TIMER.DAT 
FILEIl  =  'MODEL.DAT' 


!  Nonnal  output  file 
!  Log  file 
!  Renuis  file 
!  File  with  timer  data 
!  First  input  data  file 


*  FILED  = ' ' 

*  FILED  =  ' ' 

*  FILED  = '  ■ 

*  EOLED  =  ' ' 


!  Second  input  data  file  (not  used) 
!  Third  input  data  file  (not  used) 

!  Fourtli  input  data  file  (not  used) 

!  Fifth  input  data  file  (not  used) 
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Abbreviation 


Expianation 


Dimension 


AH(i) 

Absolute  height  of  vegetation  on  top  of  stratum  1, 
measured  from  the  plant  top 

'  AMAX 

Actual  CO2  assimilation  rate  at  light  saturation  for 
Individual  shoots 

AMDVST 

Developmental  phase  effect  on  AMX  (relative) 

AMTMP 

Daytime  temperature  effect  on  AMX  (relative) 

AMTMPT 

Table  of  AMX  as  function  of  DDTMP 

AMX 

Potential  COg  assimilation  rate  at  light  saturation  for 
shoot  tips 

ASRQ 

Assimilate  requirement  for  plant  dry  matter  production 

ATMTR 

Atmospheric  transmission  coefficient 

COSLD 

Intermediate  variable  in  calculating  solar  height 

CRRIZ 

Critical  weight  of  the  rhIzome/root  crown  system 

CVT 

Conversion  factor  of  translocated  dry  matter  into  CHgO 

DAVTMP 

Daily  average  temperature 

DAY 

Day  number  (January  1=1) 

DAYEM 

First  Julian  day  number 

DAYL 

Day  length 

□DELAY 

Integer  value  of  DELAY 

DDTMP 

Daily  average  daytime  temperature 

DEC 

Declination  of  the  sun 

DELAY 

Lag  period  chosen  to  relate  water  temperature  to  air 
temp.,  in  cases  where  water  temp,  has  not  been 
measured 

DEPTH 

Water  depth 

DLV 

Death  rate  of  leaves 

DMPC(i) 

Dry  matter  allocation  to  each  plant  layer  (relative) 

DMPCT 

Table  of  DMPC  as  function  of  water  depth  (relative) 

DSINB 

Integral  of  SINB  over  the  day 

DSINBE 

Daily  total  of  effective  solar  height 

DRT 

Death  rate  of  roots 

DSO 

Daily  extraterrestrial  radiation 

DST 

Death  rate  of  stems 

DTEFF 

Daily  effective  temperature 

DTGA 

Daily  total  gross  CO,  assimilation  of  the  plant 

g  C02.gDW-\h-' 


g  C02.gDW-\h 


g  CHaO.g  DW' 


C 


radians 


m 


g  DW.  m  ' 


g  DW.  m -^d' 


J.m’^.d’ 


g  DW.m -^d' 


C 


CO,.m^d-' 


Abbreviation 


Explanation 


Dimension 


Measured  daily  total  global  radiation 


Development  rate  as  function  of  dally  average 
temperature  sum 


Table  of  postanthesis  development  rate  as  function  of 
dally  average  temperature  sum  (used  for  calibration; 
not  read  from  MODEL.DAT) 


Table  of  preanthesis  development  rate  as  function  of 
daily  average  temperature  sum  (used  for  calibration; 
not  read  from  MODEL.DAT) 


Development  phase  of  the  plant 


Initial  light-use  efficiency  for  shoots 


Instantaneous  CO2  assimilation  rate  of  the  plant 


Instantaneous  CO2  assimilation  rate  per  depth  layer 


Leaf  dry  matter  allocation  to  each  layer  of  the  plant 
(relative) 


Table  to  read  FL  as  function  of  DVS 


Fraction  of  total  dry  matter  increase  allocated  to  leaves 


Table  to  read  FLV  as  function  of  DVS 


Diffuse  radiation  as  a  fraction  of  total  solar  radiation 


Fraction  of  total  dry  matter  increase  allocated  to  roots 


Table  to  read  FRT  as  function  of  DVS 


Fraction  of  total  dry  matter  increase  allocated  to  stems 


Table  to  read  FST  as  function  of  DVS 


Dry  matter  growth  rate  of  leaves 


Dally  total  gross  assimilation  rate  of  the  community 


Dry  matter  growth  rate  of  roots 


Dry  matter  growth  rate  of  stems 


Dry  matter  growth  rate  of  vegetation  (plant  excluding 
rhizome/root  crown  system) 


Harvesting  (0=no  harvesting,  1  ^harvesting) 


Harvesting  day  number 


Harvesting  depth  (measured  from  water  surface) 


Height  on  top  of  stratum  I  (measured  from  water 
surface) 


Selected  hour  during  the  day 


Counter  in  DO  LOOP 


Total  irradlance  absorbed  per  plant  layer 


gC02.J-' 


g  C02.m*lh-' 


g  COs.mlh' 


g  DW.m'^.d-^ 


g  CH20.m -^d' 


g  DW.m-2  d  ' 


g  DW.m*2.d-^ 


g  DW.m-2.d-' 


Abbreviation 


Explanation 


Dimension 


lABSL(i) 

Total  irradiance  absorbed  by  plant  shoots 

J.m'^.s'^ 

IDAY 

Integer  equivalent  of  variable  DAY 

d 

IRS 

Total  Irradiance  just  under  the  water  surface 

J.m’^s*^ 

IRZ(i) 

Total  irradiance  on  top  of  depth  layer  I 

J.m'^.s"* 

IWGRIZ 

Initial  weight  of  live  rhizome/root  crown  system 

g  DW.m-" 

IWLD1,2,3 

Initial  dry  matter  of  dead  leaves  cohort  1 ,2,3 

g  DW.m-* 

IWLG1,2,3 

Initial  dry  matter  of  green  (live)  leaves  cohort  1 ,2,3 

g  DW.nr'^ 

IWRIZD 

Initial  weight  of  dead  rhizome/root  crown  system 

g  DW.m-" 

IWRD1.2,3 

Initial  dry  matter  of  dead  roots  cohort  1 ,2,3 

g  DW.m*^ 

IWRG1,2,3 

Initial  dry  matter  of  green  (live)  roots  cohort  1 ,2,3 

g  DW.nr® 

IWSD1,2,3 

Initial  dry  matter  of  dead  stems  cohort  1 ,2,3 

g  DW.m-^ 

IWSG1,2,3 

Initial  dry  matter  of  green  (live)  stems  cohort  1 ,2,3 

g  DW.m-* 

K 

Plant  species  specific  light-extinction  coefficient 

m^g  DW’ 

KT 

Table  to  read  K  as  function  of  DVS 

m^g  DW',  - 

L 

Water  type  specific  light-extinction  coefficient 

m-' 

LAT 

Latitude  of  the  site 

degrees 

LT 

Table  to  read  L  as  function  of  day  number 

m'\  d 

MAINT 

Maintenance  respiration  rate  of  the  plant 

g  CHjO.m-^d-' 

MAINRT 

Maintenance  respiration  rate  of  the  rhizome/root  crown 
system 

g  CHjO.m-^'d-’ 

MAINTS 

Maintenance  respiration  rate  of  the  plant  at  reference 
temperature 

g  CHjO.m-^d' 

NGLV 

Net  growth  rate  of  leaves 

g  DW.m -^d’ 

NGRT 

Net  growth  rate  of  roots 

g  DW.nr^d-' 

NGST 

Net  growth  rate  of  stems 

g  DW.m-^d-' 

NPL 

Plant  density 

plants  .m'^ 

NUL 

Zero  (0) 

- 

PAR 

Instantaneous  flux  of  photosynthetically  active 
radiation 

J.m’^s'^ 

PARDIF 

Instantaneous  flux  of  diffuse  PAR 

J.m'^.s*^ 

PARDIR 

Instantaneous  flux  of  direct  PAR 

J.m'^.s'^ 

PI 

Ratio  of  circumference  to  diameter  of  circle 

- 

RAD 

Factor  to  convert  degrees  to  radians 

radlans.degree’^ 

(Sheet  3  of  5) 
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Abbreviation 


Explanation 


Dimension 


RC 

Reflection  coefficient  of  irradiation  at  water  surface 
(relative) 

- 

RCSHST 

Relation  coefficient  rhizome/root  crown  weight-stem 
length 

m.  g  DW’ 

RDR 

Relative  death  rate  of  leaves  (on  DW  basis) 

d' 

RDRIZ 

Relative  death  rate  of  rhizome/root  crown  system  (on 

DW  basis) 

d-' 

RDRT 

Table  to  read  RDR  as  function  of  DAVTMP 

d-'.°C 

RDS 

Relative  death  rate  of  stems  and  roots  (on  DW  basis) 

d-' 

RDST 

Table  to  read  RDS  as  function  of  DAVTMP 

d-',°C 

REDAM 

Reduction  factor  to  relate  AMX  to  pH  and  oxygen 
levels  of  the  water  as  function  of  DVS  (relative) 

- 

REDF(I) 

Reduction  factor  for  AMX  to  account  for  senescence 
plant  parts  over  vertical  axis  of  vegetation  (relative) 

- 

REDFT 

Table  to  read  factor  to  reduce  AMX  over  vertical  axis  of 
vegetation  (relative) 

- 

REMOB1,2,3 

Remobilization  rate  of  carbohydrates  cohort  1 ,2,3 

g  CHjO.m -’’.d' 

ROC 

Relative  conversion  rate  of  rhizome/root  crown  into 
plant  material 

gCH^O.gDW'.d' 

SC 

Solar  constant  corrected  for  varying  distance  sun-earth 

J.m^s*^ 

SC(i) 

Standing  crop  in  depth  layer  I 

g  DW.m'^.layer^ 

SiNB 

Sine  of  solar  elevation 

- 

SINLD 

Intermediate  variable  in  calculating  solar  declination 

- 

STEMLE 

Stem  length 

m 

TBASE 

Base  temperature  for  juvenile  plant  growth 

°C 

TEFF 

Factor  accounting  for  effect  of  daily  average  daytime 
temperature  on  maintenance  respiration 

- 

TEFFT 

Table  to  read  TEFF  as  function  of  DDTMP 

TGRIZ 

Total  live  rhizome/root  crown  system  weight  of  the 
previous  day 

g  DW.m-* 

TGW 

Total  live  plant  dry  weight  (excluding  rhizome/root 
crown  system) 

g  DW.m’’ 

TGWM 

Total  live  plant  dry  weight  measured  (field  site) 

g  DW.m-^ 

TGWMT 

Table  to  read  TGWM  as  function  of  day  number 

g  DW.m d 

TL 

Thickness  each  plant  layer 

m 

TMAX 

Daily  maximum  temperature 

X 

TMIN 

Daily  minimum  temperature 

X 

TMPSUM 

Temperature  sum  after  1  January 

X 
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Abbreviation 


Explanation 


Dimension 


TRAFAC 


TRANS1,2,3 


TREMOB 


TWLD1,2,3 


TWLG1,2,3 


TWRD1,2,3 


TWRG1,2,3 


TWSD1,2,3 


TWSG1,2,3 


YRNUM 


Translocation  factor  (relative) 


Translocation  rate  of  carbohydrates  cohort  1 ,2,3 


Total  remobilization 


Total  live  +  dead  plant  dry  weight  (excluding 
rhizome/root  crown  system) 


Total  live  rhizome/root  crown  dry  weight  of  the  current 
day 


Total  dead  leaf  dry  weight  cohort  1 ,2,3 


Total  live  leaf  dry  weight  cohort  1 ,2,3 


Total  dry  weight  of  dead  leaves  2  or  3  cohorts 


Total  dry  weight  of  live  leaves  2  or  3  cohorts 


Total  dead  root  dry  weight  cohort  1 ,2,3 


Total  live  root  dry  weight  cohort  1 ,2,3 


Total  live  +  dead  rhIzome/root  crown  system  weight 


Total  dead  rhIzome/root  crown  system  weight 


Total  dry  weight  of  dead  roots  2  or  3  cohorts 


Total  dry  weight  of  live  roots  2  or  3  cohorts 


Total  dry  weight  of  dead  stems  2  or  3  cohorts 


Total  live  stem  dry  weight  cohort  1 ,2,3 


Total  dry  weight  of  dead  stems  2  or  3  cohorts 


Total  dry  weight  of  live  stems  2  or  3  cohorts 


Dry  weight  of  leaves  (live  +  dead) 


Dry  weight  of  roots  (live  +  dead) 


Dry  weight  of  stems  (live  +  dead) 


Daily  water  temperature 


Table  to  read  WTMP  as  function  of  day  number 


Year  number  simulation  (1  -5 


g  CHgO.m'^d^ 


g  CHgO.m'" 


g  DW.m'2 


g  DW-m  ' 


g  DW.m  " 


g  DW.m'^ 


g  DW.m  " 


g  DW.m  " 


g  DW.m'^ 


g  DW.m'^ 


g  DW.m  " 


g  DW.m*^ 


g  DW.m'^ 


g  DW.m-' 


g  DW.m'^ 


g  DW-m  ' 


g  DW.m'^ 


g  DW.m  ' 


g  DW.m'^ 


g  DW.m'" 


g  DW.m'" 


Appendix  C 

Manipulation  of  Literature  Data 
Used  for  the  Model  Equations 


Photosynthesis 

Effect  of  daytime  temperature  on  photosynthesis  (AMTMP) 


To  calibrate  the  relationship  between  temperature  and  photosynthetic  activity, 
the  photosynthetic  rates  compared  with  the  photosynthetic  rate  at  35  °C  pub¬ 
lished  by  Titus  and  Adams  (1979a, b)  were  used.^ 


Table  Cl 

Relative  Photosynthetic  Activity  of  Milfoil  Shoots  in  Response  to 
Temperature  (Conditions  were  light  saturating  and  water  was  in 
equilibrium  with  atmospheric  CO,) 

Temperature, 

Relative  Photosynthetic  Rate 

0 

0.00001 

5 

0.18 

10 

0.23 

15 

0.40 

20 

0.63 

25 

0.78 

30 

0.95 

35 

1.00 

40 

0.78 

45 

0.38 

50 

0.05 

55 

0.00001 

^  References  cited  in  this  appendix  are  located  at  the  end  of  the  main  text. 
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Growth 

Assimilate  requirement  for  dry  matter  production  (ASRQ) 

The  value  of  the  conversion  factor  for  growth  of  plant  biomass,  weighted 
according  to  its  composition,  can  be  computed  in  a  simple  way  from  the  frac¬ 
tions  of  nonstructural  carbohydrates,  proteins,  fats,  cellulose,  organic  acids,  and 
minerals  (Table  C2).  This  conversion  factor  indicates  the  amount  of  glucose 
consumed  to  produce  each  g  of  plant  biomass  (g  CH2O  g  DW^).  This  method  has 
been  employed  to  calculate  assimilate  requirement  of  milfoil  shoots  for  biomass 
production. 


Table  C2 

Estimated  Chemical  Composition  of  Milfoil  Shoots  (this  study)  and 
Typical  Conversion  Efficiencies  for  Agricultural  Crops  Showing 

How  Much  Glucose  is  Used  for  the  Synthesis  of  Each  Organic 
Matter  Component  (Penning  de  Vries  and  Van  Laar  1982b) 

Component 

Contribution  to  Biomass 
percent 

Conversion  Factor 
g  CHjO  g  DW’ 

Nonstructural  carbohydrates 

14 

1.242 

Proteins 

17 

1.704 

Fats 

8 

3.106 

Cellulose 

33 

2.174 

Organic  acids 

11.2 

0.929 

Minerals 

16.8 

0.050 

Milfoil  shoot 

100 

1.539 

Note:  As  the  conversion  factor  for  cellulose  was  not  known,  that  for  lignin  has  been  used. 

Site-Specific  Environmentai  Conditions 

pH,  alkalinity,  and  trophic  state 

pH,  alkalinity,  and  trophic  state  are  important  factors  influencing  primary 
production  in  aquatic  systems.  pH  and  alkalinity  determine  carbon  availability 
for  photosynthesis,  and  trophic  state  gives  an  indication  of  algal  production  and 
consequent  light  attenuation  within  the  water  column.  The  model  is  calibrated 
for  dissolved  inorganic  carbon  concentrations  1. 1-1.8  mmol  (alkalinity  Lake 
Wingra  1. 1-1.8  mmol;  Lee  and  Kluesener  1972).  pH  affecting  potential  photo¬ 
synthetic  rate  at  light  saturation  through  REDAM  can  be  modified  by  the  user. 


C2 
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The  model  is  calibrated  for  a  light-extinction  coefficient  range  of  the  water  of 
1.15  -  2.0  m  ‘  (Lee  and  Kluesener  1972);  the  value  of  this  parameter  (L)  can  be 
modified  by  the  user. 


Water  temperature 

The  temperature  has  been  measured  in  the  surface  water  of  Lake  Wingra  at 
several  points  in  time  in  1970.‘  For  Days  1  and  365,  the  same  temperatures  as 
those  measured  on  the  nearest  dates  in  Lake  Wingra,  Wisconsin,  have  been 
taken. 


Table  C3 

Seasonally  Measured  Daytime  Temperatures  in  the  Surface  Water 
of  Lake  Wingra,  Wisconsin,  during  1970 

Day,  number 

Temperature,  ""C 

Day,  number 

Temperature,  ""C 

1 

3.5 

216 

25.3 

62 

3.5 

223 

25.3 

69 

4.0 

230 

24.4 

76 

5.3 

237 

23.2 

84 

6.3 

244 

22.5 

90 

6.5 

246 

22.9 

97 

5.7 

251 

23.0 

98 

7.0 

258 

16.5 

104 

6.7 

265 

20.0 

111 

6.9 

272 

15.8 

118 

15.2 

278 

15.1 

125 

15.3 

286 

14.3 

132 

17.6 

293 

11.8 

139 

17.0 

300 

12.8 

146 

19.1 

307 

8.2 

153 

19.1 

321 

4.1 

160 

22.7 

328 

0.3 

167 

23.9 

335 

1.7 

174 

22.7 

342 

0.9 

181 

24.8 

349 

0.1 

188 

23.5 

355 

1.2 

195 

26.8 

363 

1.6 

202 

22.4 

365 

1.6 

209 

26.7 

‘  Personal  Communication,  1995,  J.  E.  Titus,  Univeristy  of  Binghamton,  New  York. 
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downward  transport  of  soluble  carbohydrates  occurs  after  anthesis  of  each  cohort,  replenishing  the  rhizome/root 
crown  system.  In  a  tropical  climate,  a  third  plant  cohort  is  active. 

MILFO  simulated  the  dynamics  of  plant  and  rhizome/root  crown  biomass  at  Lake  Wingra  well  over  a  period  of 
1  to  5  years.  It  has  been  used  to  calculate  plant  and  rhizome/root  crown  biomass  for  the  same  latitude  in  a  different 
year  and  for  other  latitudes  in  temperature  (Alabama,  USA)  and  tropical  (India)  areas,  where  it  simulated  biomass 
ranges  similar  to  those  measured  in  the  field. 

Sensitivity  analysis  showed  that  maximum  plant  biomass  of  a  Eurasian  watermilfoil  community  is  most  sensitive 
to  a  change  in  photosynthetic  activity  at  light  saturation  and  very  sensitive  to  a  change  in  light-use  efficiency,  and 
that  end-of-year  rhizome/root  crown  biomass  was  often  more  sensitive  than  maximum  plant  biomass.  The  latter 
illustrates  the  utmost  importance  of  the  rhizome/root  crown  system  for  local  survival  and  biomass  production  in 
milfoil. 

Environmental  factor  analysis  indicated  that  changes  in  climate  can  greatly  affect  simulated  end-of-year 
rhizome/root  crown  biomass.  Maximum  plant  biomass  proved  far  more  sensitive  to  changes  in  water  transparency 
than  to  changes  in  water  depth. 

The  model  can  be  used  as  a  tool  to  predict  the  dynamics  of  a  Eurasian  watermilfoil  community  over  1-  to  5-year 
periods.  Running  the  model  with  different  parameter  values  specific  for  any  particular  site  and/or  treatment,  for 
example,  biomass  removal  to  a  certain  water  depth,  helps  in  gaining  insight  into  the  predominant  mechanisms 
regulating  submersed  plant  dynamics. 


